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PREFACE 


The work presented In this final report was In response to the following 
three topics as suggested In the contract's scope of work. 

1) . Investigation of possible multivariate extensions of existing 

univariate distributions which have been used for modeling 

meteorological phenomenon. 

2) . Development of Goodness-of-fl t tests. In particular for non- 

Gaussian distributions. 

3) . Investigation of the effect of correlated observations on 

statistical inference 

Reports 1-4 are concerned with some aspects of topic #1. Report 1 contains 
an estimation procedure for several discrete multivariate distributions. 
Report 2 contains a procedure for computing cloud cover frequencies in the 
bivariate case. This procedure can be used to compute probabilities for 
cloud frequencies fcr either two geographical locations or for the same 
location at different times. Report 3 contains the procedure and correspond- 
ing computer code for calculating conditional bivariate normal parameters. 
This report was requested by the COR. Report 4 contains a procedure for 
transforming multivariate non-Gaussian distributions into a nearly Gaussian 
distribution. 

Reports 5 and 6 are concerned with topic #2. Report 5 contains a 
goodness-of-fit test for the extreme value distribution which is used in many 
meteorological applications. Report 6 contains a goodness-of-fit test for 
several continuous distributions. 
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Report 7 Is concerned with the problem given In topic #3. In this 
report, the effect of autocorrelated observations on confidence regions Is 
Investigated. 

Report 8 contains a computer code for generating both random and non- 
random observations for specified distributions. This program was used to 
generate the samples for the Monte Carlo simulation needed In the other 
reports. 
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Summary 



Procedures for estimating the parameters of three discrete 
multivariate distributions, the Multinomial, Negative Multinomial, 
and the multivariate Poisson distribution, are given along with 
approximate variances for the parameter estimates. 


I . INTRODUCTION 


This paper is concerned with the problems associated with the 
estimation of parameters for three discrete multivariate distributions, 
the multinomial, negative multinomial, and the multivariate Poisson, 
which are the multivariate extensions of three common univariate 
discrete distributions, the binomial, the negative binomial, and the 
Poisson distribution. The distributions are introduced in Section 2. 

A detailed explanation of the estimation procedures along with 
approximate bounds for the variances of the estimates are given in 
Section 3* An example is presented in Section 4 which is intended 
to demonstrate the use of the estimation procedures. A listing 
and card input description of the computor program is given in the 
Appendix . 
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II. DISTRIBUTIONS 


Johnson and Kotz ( 1969 , Ch. 11) provides a detailed 
discussion of the functions described below. 

2.1 Multinomial Distribution 


The simplest of the three distributions both in 
structure and theory is the multinomial distribution. Let 
E^^,...!^ be possible events which can occur from a series 
of independent trials. If E. has probability P. of occuring 

U u 

and nj is the number of times E . occurs in the N trials where 
o 0 

k 

£ n . = N , then the joint distribution of the random 

3=1 3 

variables n., ,n 0 , . . . ,n v is the multinomial distribution with 
parameters N,P^,P£, . . . ,P^. The distribution is defined by 

^ n i k 

PCn^n^ . . . ,n k ) = N! n (P^ 3 /n..!) (91^; z . n ^ =N )* 

j = 1 0 “ 1 


( 1 ) 


2.2 Negative Multinomial 

Just as the multinomial distribution is a natural 

extension of the binomial distribution, the multivariate 

negative binomial distribution is a natural extension of the 

negative binomial distribution. Hence, the probability generating 

function for the multivariate negative binomial is defined by 

k M 

(Q - i P 4 t .)" N 
i=l 1 x 


( 2 ) 
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with P* > 0 for all i-l,...,k; N > 0, and Q - i P. « 1. 

i-1 1 

Prom formula (2) we have the following distribution function 

k 

r(N+ ? n i> k n. 

P(n*^ >np , • • • ^n^ ) * k Q n (P^/Q) ^ (3) 

( n n. !) (N) ‘ i = 1 

i-1 1 

where n i > 0, i-1, . . . ,k. 

This is called the negative multinomial (or multi- 
variate negative binomial) distribution with parameters 
N ’^l’*2’ * ’ * ’^k 1 where N is a non-negative integer. A special 
form of this distribution is a compound Poisson distribution 
which can be further simplified to a bivariate form as 
described by Bates and Neyman (1952). 


2.5 Multivariate Poisson 


Consider a sequence of k variables x^^,... ,x k 
such that each one is a combination of two independent uni- 
variate Poisson variables where one of the Poisson variables 
is present in all k variables. That is, 

x i = u+v i» x 2 = u+v 2 ’ • • • » x k = u+v i and u ’ v l’ v 2’ * • * ’ v k 
are independent univariate Poisson variables with expected 
values 5 » e ii e 2’-“» e k respectively. The Joint distribution 


of , x^ i • • • i X k ^ 

■2(x^ * • • • ,x k ^ a ®xp(— £— ‘ 


m 
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-V 

x 2 -d 

_®2 

Txp-j)! 


d=o 


x x -d 

i 0 x 
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dT ~7xp3TT 
x k -d 

e k 

Tv-'d ) ! 


• • • 
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where m - min(x 1 ,x 2 , . . .x k ) . This is called the multivariate 
Poisson distribution with parameters C O 2 * • • • » 6 k* 


III. ESTIMATION 

In section 3.1 the techniques used to estimate the 
parameters of the three distributions are described. The 
subsequent section is concerned with the variances of the 
estimates for the multinomial and negative multinomial dis- 
tributions. A computer program was written to perform the 
needed computations. 

3.1 Parameter Estimation 

The maximum likelihood estimates of * P 2 •> * * for 
multinomial distribution are the relative frequencies 




n d/N 


(j=l, ... ,k) 


(5) 


where n. is the observed frequency of E. given N independent 
trials. 

The method of moments is the most convenient approach 
for estimating the parameters of the negative multinomial 
distribution. The moment generating function of a k variate 
negative multinomial distribution is 

k -N 

m (t^ j j t k ) - (Q — r P^c ) • 

i =1 


Thus we obtain the following moments 

3m(tj t . . . ) 


E(n d ) 




= NP^ for J=1 , . . . ,k 


t=0 
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E(n.nj) 


iV * 


3 ) 

«r - 


t«o 


N(N+l)P i Pj 


giving 


and 


• I,2p i ? i j +NP i I> d 

E(n. )E(n j ) 
« E(n i ) E(n^) + -j p -* 3 - 

E(n i )E (n^) 

N " tTn in "j-)“E(n i )ECn d T 


(6) 


Pj« E(n^)/N. (7) 

Equating raw estimates to moments to obtain an estimate for N, 
we have 


N « 


n i n 1 


n i n j “ n i n j 


and P. = n^/N for i, j=l,...,k and iYj where 

O J 


n n 

n^ « -ill — i- ; n^n^ = 1= - 1 n 1 — - given n observations. 

The accompanying computer program utilizes this method of 
moments in two ways. There are k(k-?.)/2 possible estimates of 
N by this method where k is the number of parameters. Similarly 
there are k(k-l)/2 possible values of n.n . as well as n.n .. The 

1 J ^ U 

program first average j the k(k-l)/2 values of n^n.. and n^n^ and 
then outputs an estimate of II based on these averages. The 

second approach calculates the k(k-l)/2 estimates of N and prints 

out the average estimate of N. The parameters P^,i=l,...,k is 

also estimates twice corresponding to the two estimates of N. 
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The method of moments is also used in estimating 
the parameters of a multivariate Poisson. The moment gener- 
ating function is given by 

T k k 1 

m(t,,...,t,.) - exp -c(l-exp( I t.)) - z e.(l-e ) (f 
A K L i-1 1 i-1 J 

It follows that 




ECxj) + e. 


^m(t 1 ,... ,t k ) 




Therefore 


- (e+e i )(5 + o J ) + c 

* E(x i ) E(xj) + 5 


5 ■ E [*i x a] ' E [ x i] E [ x a] • 


Substituting raw estimates for expected values we have 


n 

* x t . 

t-¥j- ; i x a where *r i= V i •- ¥k 


1 X, x v . 

i-1 Kl . 


(ID 


Since 8^ * E(x^)-c, a method of moments estimate for is 

^ A 

(Xj-c). Again the accompanying computer program uses two 
approaches to estimate C via the method of moments. First 
the program averages all possible values for x- x . and x.x, and 

J J 1 J 

estimates $ based on these two averages. Next the program 
averages the k(k-l)/2 possible estimates of C and outputs 



this average as a workable estimate of £ . The parameters 


of e^, i=l,...,k are estimated twice to correspond to the 
two estimates considered for 5 • 

3.2 Variances of Parameter Estimates 


The exact variance of the estimates for the multinomial 
t rameters can be easily derived. Consider 

a 

var (P.) = var (n./N). 
tj <J 

= E(n./N) 2 - {E(n ./N)} 2 
J J 

4 - (N 2 P 2 + NP.q.) - P. 2 

t J 


N' 


p i q i 

-iP 


■ J it J 


( 12 ) 


hence an approximate variance for P. is P.(l-P_. )/N. 

u J J 

In order to place approximate bounds on the variances 
of the negative multinomial parameter estimates, consider 
Fisher's Information Matrix for the maximum likelihood para- 
meter estimates which is defined as 

-1 


V(a-^ , • * • , ) — (E 


3 a . 9 a . 
- ^ J 


) 


(13) 


where a- and a. are parameters and L is the likelihood function. 
-L J 

Kendall and Stuart have shown that this matrix is the asymp- 
totic variance-covariance matrix for the maximum likelihood 
parameter estimates. From equation (3), we have the following 


n k 

n r(N+z n. .) 

L = Jr 1 1J 

3-1 Li ni i ! )( r W) n 
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-Nn 


n (Pj^/Q)^' 
i = l 


.i n « 


(i4) 
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In L 


n n k 

= s In r(N+Sj) - In (n n n. .,!) ~ n lnr(N) (15) 

0=1 ,1=1 i = l 0 


k n 


-Nn 


In Q + E (Z n. . )(ln P. - In Q) 
i=l cl 10 1 


k 


n 


where S. = E n. S. = E n. . , n is the number of samples 

' n i=l 10 1 j = l J 

taken and n. . is the number of times E. is satisfied on the 


0 th sample. 
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3 In L 
3 N 


n 2(3. )-l 
= E E J 
0=1 k =0 


N+K 


- n In Q 


(16) 


3 In L 
3*-N 


n E(S . )-l 

= E E l ' 

0=1 k=0 


A 

— =1_5 = Z - (N+k-l)~ 2 E(F.) (17) 

(N+k)^ k=l J 


where F. is the number of S.'s greater than or equal to 0 , 

0 0 


3 2 ln L 

~mr~ 


-n 


k 

1 + Z P. 
i = l 1 


'or i=l , . . . ,k 


(18) 


3 On L 
' 3 P"." ■ 


-Nn - ) E(s l , 

+ — ■ t; — 

1 + z p. 1 

0=1 0 


(19) 


3 2 In L 


K 1 

Nn 4 z E(S p ) 

1=1 1 


(1 + * F<> ) 

t-1 


2 


E(S 1 . ) 

u 

P.^ 

0 


( 20 ) 
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Ur, inf; 

valuer, 

bounds 


Nil i » KU- . 1 'i 
3' In L l -I 1 

k . 

01 ( l -*■ e r, V 

l --i 1 

the two nets of erst; i mute:; for 

for I-J(S .), ECS . 1 ), and E(F.) 
• * d 0 

for VCN,^,... ,f k ). 


• for i / ,i (SI ^ 

N, r x , - • • ,E k and numerical 
wo can obtain approximate 


IV. AN EXAMPLE 


Negative multinomial data were obtained from Arbous and 
Kerrich (1951»P* **24) to illustrate the output from the computer 
program. The results are found in Table 1. Notice that in the 
binomial case both estimates of N are the same since there are 
only two variables. For this same reason, only one Fisher’s 
infoi’ination matrix is produced. If more than two variables 
were considered, we would have obtained two different estimates 
for N and the information matrix. From the two distinct varian- 
ces obtained from these matrices one could obtain the boundary 
points of the internal about the variance of the parameter 


estimates . 
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TABLE 1. 


THE MOMENT ESTIMATE OF N OBTAINED BY AVERAGING 

THW RAW MOMENTS FIRST IS 3-350 

THE CORRESPONDING PROBABILITIES ASSOCIATED WITH 

THE RESPECTIVE VARIABLES ARE 0.295 

0-385 

THE MOMENT ESTIMATE OF N OBTAINED BY AVERAGING 

ALL POSSIBLE MOMENT ESTIMATES 3-350 

THE CORRESPONDING PROBABILITIES ASSOCIATED WITH 

THE RESPECTIVE VARIABLES ARE 0.295 


0-385 

FISHER'S INFORMATION MATRIX USING THE MINIMUM 
ESTIMATE OF N 

1.207 

- 0.108 0.010 
-0.140 0.012 0.017 
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APPENDIX 


Cols. 

3 


1-3 

4-7 

7-77 

1-70 


CARD INPUT DESCRIPTION 
Card 1 

1 if a multivariate poisson distribution is to be analyzed 

2 if a multinomial distribution is to be analyzed 

Any other number in this column indicates that the 
negative multinomial distribution is to be analyzed. 

FOR THE MULTIVARIATE POISSON AND 
NEGATIVE MULTINOMIAL DISTRIBUTIONS 

Card 2 

contains the number of variables 
contains the number of observations 

contains 7 pieces of data in consecutive 10-column spaces 

Card 3 + 

contains 7 pieces of data in consecutive 10-column spaces 


FOR THE MULTINOMIAL DISTRIBUTION 
Card 2 

1-3 contains the number of events 

4-74 contains 7 pieces of data in consecutive 10-column spaces 

Card 3 + 

1-70 contains 7 pieces of data in consecutive 10-column spaces 


APPENDIX 

IMPLICIT REaL*9 U-H/D-ZI 

Rc Al*8 Ell 2 l#Efc l» tPf li I#NJ« 35-» 12 ANl .4* )#crt( 

cL)/TI12)/FI3Sv>>/ INP(:3/13I/S<35CI/SPI12I/RMHI> 

28 RcAOI5/ 1/eNOU' I Irt 

IF I IH.EQ.2I GO TU 30 . , v . , , Ui 

READ 1 5/ a I K/M/KNJI l/JI/J-i/M/I-l/NI 

1 F0RMATt2I3/l7F:;»..l> 

Ki-K^l 

KK-K-1 
00 9 i-i/K 
El 11-000 
9 SPI I l-ODC 
00 *0 I -1/ KK 
00 10 J -1# KK 

10 EP 1 1 / J 1-000 

u Hitei! 1 '" 

DO 14 I -1/Kl 
DO 14 J-i/Ki 

14 inf i i# ji-odc 
00 2 I-l/K 
DO 3 J-l/M 

3 El II -El II4NJI J» 1 1 

2 El 1 1-El Il/M 
DO 4 J-l/KK 
JJ- <!♦! 

00 4 L-JJ/K 
Ll-L-i 
DO 5 l -1 / M 

5 EPI J/Lll-EPI J/LU+NJI 1/ J) *NJU#l I 
EE I J/ltl-EtJI*elll 

4 EPI J/LLI-EPI J*LU/1 

51 - LOO 

52- GDC 

00 b I -1 /KK 
00 b J • I /KK 
Si-Sl+EEI I / J ) 
b S2-S2*EP( 1/ J I 
G-K* IK-1 I/2DC 

CALCUlAn ON^F* THE NE oi ° MU LT I N'JM I AL PAf\ A ME TF k $ BEGINS HERE 
HN-S1/ IS2-S1 ) 

IF ( HN.LE.ODC I GC TO 60 

11 WRlTEIb/0) HN 
WRITE (o,26) 

DO 7 I-l/K 

pi n-Ei ii/hn 

7 WRITEI6/27I P( I I 
SUM-000 
DO 13 I-l/KK 
DO 13 J-I/KK 

AN I I / J I > E E I I / J I / I : P I I » J I ~ ^ * I/JII 
IF lANII/JI.LE.l'Ow) GD TO su 
13 SUM«SUM+AN( 1/ J I 
SUM-SUM/G 
WRITE lb/241 SU* 

WRITE (b/^6) 

00 15 I-l/K 
PI I )-E( ll / SUM 

15 WRITE! b/27 I P( I I 
0- AN 1 1 / 1 1 
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33# 33# 32 
‘ MULTINOMIAL 


tlSTRlUUTION 


Ae*ANt i#l) 

DO 5 L I •£# KK 
DO 31 J«I#KK 
IF (AN(1#J)-C> _ 

N>C IN TH£ NEGATIVE 

32 O a AN ( I # J ) 

33 IF I AN! I# J )*AE ) 3<»#3a#31 

34 Ae-AN( I# J) 

3 i CONTINUE 

D IS NOW THE NAX ESTIMATE OF N ANO Afc IS THE MIN 
Ob 62 I *1# K 
00 o2 J-1#M 
SP( I I • SP ( I > ♦N J ( J# 1 1 
DO 63 I«1#N 
DO 63 J»l# K 
S( I ) B S ( I 1 ♦ N J 4 I # J I 

DO 64 I-l#MM#i 
II* l *2 

I) 75#7s#7b 

DD-sTi ) 

S( II'SIIM 
S(I*il-DC 
DO 64 J« II# H 
IF (SU)-SUII 3i#o5»34 

oo»sm 


62 


63 


75 


76 


65 


64 


6d 

67 

66 


78 

36 

35 


36 


37 


38 




S( ID-00 
S( J)*CE 
CONTINUE 
L-S(H) 

00 66 J»1#L 

Du 67 I»i#M 

IF (SC I)-J> 67# fc 3# bo 

F( JJ-H-D100 

GO TO 66 

CONTINUE 

CONTINUE 

SH-bDC 

00 78 I"l# M 

SH*5H+S( 1 1 

DO 35 I • i» L ..... 

INF( 1#1 )*INF(i# 1 ) ♦( i Cb/ (O + I- .tC ) **2 ) *F ( I) 
SP I *CD J 
DO 36 I *1# K 


PCJI-£|IJ/D 


SPI + PU ) 

00 37 !»2#Ki 
H«M 

INF ( 1# I) «H/( IDl +SP I ) 

DO 38 J«2» KI 
00 38 I • J# Kl 

INF ( J# 1 1 ■- ( D-M + SH) / ( i 00 ♦SP I ) **2 
IF (I.fcQ.JI INF( I# J»- IN F ( I. J > ♦ S P f J- 
IF (D.EO.AE) WRITE (o»*A) 

IF (O.NE.AEI WRITE (b»*3l 
CALL ARRAY <2#K1»RN# INF ) 

Call sinv <rm#k:#.05# ieri 

CALL ARRAY(l#Kl#RM# INF) 


I/PIMIM2 


'V 

* ^my 


) 
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DO 23 I-l.K 
li *1 ♦! 

00 23 J-Ii.Kl 

23 INF I J, 1 1 - I NF ( 1/ J) 

00 4l J* 1, K 1 

41 WRITE (6. 2) ( 1 NF I I. J). !■*# J) 

42 FORMAT (9F14.8/4F1**. 8) M 

43 FORMAT I • CF I SHERS INFORMATION MATRIX USING THE MaXIMJM ESTlMATe CF 
CN* I 

IF (O.EO.AE) GO TO 28 
O a AE 

GO TO 34 

44 FORMAT('CFtSHfcR'*S INFORMATION MATRIX USING THE MINIMUM ESTIMATE 
COF N * I 

CALCULATIONS OF MULTIVARIATE POISSON PaRaMETFxS 8EGIN HERE 

16 B‘IH:h!68c. od n ,o 

SUM-ODO 

WRITE (6.18) PE 
WRITE (6.21) 

00 20 I • 1. K 
T( 1 1 -E ( 1 1-PE 

20 WR IT E( 6. 27 ) Tt I ) 

DO 17 I-l.KK 

DO J-I.KK 

EH( I. J)-EP(1. J )-£E( I. J) 

IF (EH( I# JI.lE.:-Cw ) GO TO V. 

17 SUM.SUMaEH( I. J ) 

SUN-SUM/G 

WRITE (6.19) SUM 
WRITE (6.21) 

00 22 1*1. K 
T( I ) *E ( I )- SUM 
22 WRITE (6.27) T< I ) 

27 F0AMAT(31X.F:-.a) r , , ra , 

19 FORMAT ( • O' » ' THE MOMENT ESTIMATE OF THE POISSON PARAMETER . FOk U C8. 
CAINEO BY AVERAGING ALL*/* POSSIBLE MOMENT :STlMATcS IS ' » iX. F . 4. 6 ) 

8 FORMAT ( *0* . • THE MOMENT ESTIMATE OF A OBTAINED 8Y AVERAGING ThF man 
C MOMENTS FIRST l S ' / 3 i X, FI * . o > ... tcco 

26 FORMAT ( ' THE CORKt S PONO ING PKC BAB I L l T 1 6S ASSOCIATED WITH THc RESP; 
CCTIVE VARIABLES ARE*) „ ir 

24 FORMAT ( 'O' , ' THE MOMENT ESTIMATE OF A OBTAINED BY AVERAGING ALL POS 
CSIBLE MOMENT ESTIMATES'/* OF N IS' »2*x.Fl*.d) 

18 FORMAT! 'O'. 'THE MOMENT ESTIMATE OF THE PCISSUN PAR AMET tR FOR U H8T 

CAINEO BY AVERAGING THE RAW MOMENTS FIRST IS'. Fit. 8) , ,, 

21 FORMAT (' THE CORRES PONDING ESTIMaTi CF THE POISSON PaRAMETcR FOR 
c the RESPECTIVE v VARIABLES aFE'» 

GO TO 26 

CALCULATION OF MULTINOMIAL PARAMETERS BEGIN HE&E 
30 READ (5,49) K, ( T ( I ) , I *1 ,K ) 

49 FORMAT (13, ( 7F1 w ••■ ) ) 

WRITE (6,52) 

SUM* ODO 
DO 50 1-1, K 

50 SUM*SUM*T( I ) 

DO 51.1-1, K 
P( I > -T < 11/ SUM 

SP( I )-?(!)♦( IDO-Pt I ) )/SUM 

51 WR I T E ( 6, 48 ) P( I),SP( I ) 

|2 pgRMA| CJ^PRQ|AS^LITIES^>1CX,*APPRCXIHAT= VARIANCES') 

GO TO 28 . 

•’ ISiI,ilS^* l NON-NiS»TI»E ,4«*HETER H»S »EEN ESTI-UTSO »S N?GATIV£M 

stop 

END 
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SUBROUTINE ARP AY IMOD E» N #RM# INF ) 
IMPLIC IT REU*B (A-H#0-2) 

Rt AL*8 RHI9: >, INF(13#.3I 

ip inode-:) i U'j 


IJ-0 

OU lie K-i> 
00 llo l-i# 
J J - 1 3 ♦! 
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SUBROUTINE SlNVU#N»cPS#URI 
IMPLICIT REAL* o U-t<>>2) 
DIME NS (ON A(4. ) 

Call mfsci a*n# lPS* if* » 
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END 
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A PROCEDURE TO PREDICT CLOUD COVER 
FREQUENCIES IN THE BIVARIATE CASE 


Summary 


The purpose of this report is to present a procedure 
for approximating cloud cover probabilities for two different 
locations or for the same location at different times. In 
addition a monte carlo procedure Is presented for integrating 
the bivariate normal distribution. This program is used for 
computing the approximate probabilities. 

If one assumes that the density function for the bivariate 
cloud cover model Is approximately bell-shaped, then it is shown 
that the desired conditional probablities can be approximated 
using the bivariate normal distribution. Examples illustrating 
the feasibility of this procedure are included. However, if 
the bivariate density for the cloud cover model is highly J or 
U shaped this procedure provides results which are less than 
satisfaefory . Examples illustrating this situation are also 
included . 


I. INTRODUCTION 

The purpose of this report is to present a procedure 
for estimating Joint probabilities for the degree of cloud 
cover over two regions or one region at subsequent time intervals. 

Falls (197*0 demonstrated that the beta distribution 
adequately describes th' 1 variation in the amounts of cloud cover. 
This conclusion was based upon analysing cloud cover data frbm 
diverse locations, for different times of the year and for 
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different times of the day. Thus, we may expect that the 
multivariate beta distribution, sometimes called the Dirchlet. 
distribution would be a natural extension for describing the 
bivariate case. However, a theoretical requirement of the Dirchet 
distribution is that the variables be negatively corrected, and 
this constraint seems to intuitively disagree with the actual 
situations. Consequently, a different approach was required, one 
allowing for both positive and negative correlations. 

Peizer and Pratt (1968) provide a possible approach, that 
of using the norial distribution for approximating tail proba- 
bilities in the seta distribution. Thus, if one assumes that 
the correlation between the two sites is structurally related 
to the co -relation present in the bivariate normal distribution, 
one may be able to extend the work of Peizer and Pratt to the 
multivariate vetting, that of approximating Joint probabilities 
using the bivariate normal distribution (BVN). This approxi- 
mation would apprar to work adequately for those cases where the 
univariate normal approximation gives satisfactory approximations 
to the beta distribution. 

This report consists of three main sections. The first 
section describes a program for integrating the BVN over rect- 
angular regions. This section is basically self contained, and 
.It provides the user the needed explanation for integrating the 
BVN, The second section illustrates how this procedure is used 
in approximating the bivariate cloud cover model. Applications 
and examples of this procedure are presented in section 3* The 
program documentation and listings are presented in the Appendix. 



II. BVN PROGRAM 

A procedure was required for integrating the bivariate 
normal distribution over a specified region. The BVN program 
provides an approximation to the above integral. This section 
consists of three subsections , 1) introduction to the monte 
carlo theory, 2) application of this theory to the BVN distri- 
bution, 3) examples. 

2.1 General Monte Carlo Technique 

An excellent summary on the general principles of monte 
carlo theory can be found in Newman and Odell (1971). The 
following is a discussion of this method as related to double 
integration. 

Let x«(x^,X 2 ) denote an arbitrary two dimensional vector 
and f(x) a real valued function of x* Consider the integral 

m m 

• “ J J f(x)e(x) dx i d *2 (2.1) 

where g(x) denotes a probability density function on the plane. 
The integral (2.1) is the expected value of f(x) and can be 
estimated by 

. 1 N 

e * H l r(x ) 

i-i ~ 

where x^ , i*=l,...,N are random samples from the pdf g(x). The 
variance of 0, is given by 
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var (•) ■ J var I. L m dxjdx 2 

which can be estimated by 

“ *) 2 . 

n-i iml i 

The estimated standard error is given by, e = s//n. 

The following describes a procedure for reducing the 
magnitude of the var ( e ). Suppose that there exist n a function 

O 

h(x) on R' (two dimensional reals) which approximates f(x) on 

p 

R and suppose that 



h(x) g(x)dx x dx 2 


is known. Then 





(f(x) - h(x))g(x)dx 1 dx 2 . 



The variance of f(x) - h(x), is given by 

var (f(x)-h(x) ) - var (f (x))+var(h(x))-2 cov(f (x),h(x)). 

If var (h(x)) < 2 cov (f(x), h(x)), we have that 
var (f(:c) - h(x)) < var (f(x)). 

Note that if (f-h) and h are positively correlated then var (f-h) 
is less than var (f). This is true since 

var(f ) » var [ h+f-h ] 

- var(h) + var(f-h)^ cov(h,f-h) 

Thus we have 

var (f-h) ■ var(f) - var(h) - 2 cov (h,f-h). 
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Assume the correlation of (f-h) and h is positive. Hence, 

var(f-h) < var(f) - var(h) 
which implies that 

var (f-h) < var(f ). 

Therefore the larger the correlation of (f-h) and h, the 
greater the reduction of the variance by removal of the regular 
part h(x). 

2.2 Program Explanation 

The object is to integrate 

8=/ 1 / 2 f(x|jj, E)dx 1 dx 2 . 
a l a 2 

o P°x °2 

where ji'= (v^iig); * = a 2 )and 

f(xj ji, E ) = BVN distribution - 


2*0x0. 


2 1/2 
(1-P) 1/2 


exp{- 2(l-p^) 


r ,UzlLl 
L ( o 

l 


, 2 
) - 


(Xx-Qx )(x 2 -o 2 ) 


( 


x^ )2] 

02 


( 2 . 2 ) 


In formula (2.1) we define g(x) = 1 , i.e. 

^l“ a i ^ Cb 2 ~a 2 ) 

g(x) represents a bivariate uniform distribution, and 
evaluate the integral 


dx-i dx. 


/I , u O 

0 = / / f(£l il, E) T b ~ -Tl ' (b~-a V 

a x a 2 ^D 1 -a 1 AD 2 a 2 ; 


( 2 . 3 ) 
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It follows that 

• - PCbj-a^Cbg-ag) 

and the estimate 

5*5 (b^apCbg-^) 

where 

1 N 

0 * ft 1 f (iil iL» O 

i=1 -i 

when ^ is a random vector from the pdf 

8<s> 'Av ai x b 2 -a 2 ; i V x J ib j ’ r 1 - 2 

\ 0 ; Otherwise 

Since g(x) is the product of two independent uniform distributions, 

a random vector is generated using the equations x^=a .+u.(h.-a .), 

j=l,2 v/hcre u. is distributed uniform over the interval (0,1). 

J 

In the BVI1 program the regular part h(x) is defined to he 
all the terms up to the cocffeciont 1/8! in the two dimensional 
Taylor's expansion ( Folks 1989, p. 260 ). The two dimensional 
Taylor's expansion Liie point ( a ^ , a ^ ) is given by 

j IT 

f(x 1 ,x 2 ) = f(a 1 ,a 2 ) + (x 1 -a 1 ) 

O 

+ < x 2~ a 2> \^q (a l’ a ? ) + 7T ^ (x l _n l )2 <a l’ a 2 ) 

+ 2(x 1 -a 1 )(x 2 -a 2 )|;^ + (x 2 -a 2 ) j-^s (n^ 2 ) ~ 

+ 7^-f (ni. a 2 ) + 3 (x 1 -o 1 ) 2 (x 2 -o 2 )^/— ( ttl ,a 2 ) 

L ^i l 2 

+ 3(x 1 -a- L )(x 2 -a 2 ) 2 — 2^ x 2“ n 2^ ■ — £ ^ a l’ a 2^1 + 

^1^2 ^ v ^ — * 
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Hence it was necessary to find all pnrtiuls (up to 8th order) 
of the BVN distribution function, f(x^,x 2 ). 

Let (a^,a 2 ) = (w^, the mean vector of the BVN 
distribution. Then equation (2.4) becomes 


ax/ w l ,M 2 )=f(x l’ x 2 

X 9 


>r— 

- 2 ( 1 -: 




2 

T2 


( x l-wi) 


- -2£_ 
° 1°2 


(x 2 -w 2 )} 


«= 0 

x l* W 1 
x 2= u 2 


a 2 f 


iX, 


2 = 


-1 


2wo l °2^ 1_p ' 


17 ? 


a 2 f 


Sx,3x 0 2 ?/, ?n-V? 

1 ' 2*°-. o-, (1-p ) 


1 


O 

Hi; 

3x^ £ 


-1 


2 io^o 2 ' (1-p ) 


2 N 5/2 


a^f 
— 1 

dXn 


2va 1 ^a 2 (l-p 2 )' ?/2 


a 4 f 


a 


1 


. - ? P - 


2.. 1 \ 2 2 ( 1 -p 2 ) 5/2 


2U 


3*f 

i 


£ ?p ? fl 




3 4 f 


3 


-=2fi. 


a^Vci -’ 2 ) 572 


3 


3 X- 


2 * 0 1°2 5(1 " P 2)5/2 


3 6 f 


*X, 


J=ll 


2»Oi7o 2 ( ]L _ p 2 ) 7/2 


3 6 f 


3x^^3X2 


15 p 


P __ 6„ 2 n 2x7/2 
2 °1 a 2 ( 1_p ' 


3 6 f 

^7 


-3-l3>' 


2 *Oi5o 2 3(i-p 2 ) 7A 


3 6 1 


5 5 

3x^X2^ 


9p+6p 


4 4, 2x7/2 

2 * 1 °2 ( 1_p ) 


3 6 f 

: — 2 — 5 

SX-J^ 3X 2 




x 5, ?.7/2 

2 WO l^°2^( 1 ~ P ^) 


3 


> Xj» x 2 - 


15 » 


2.. 2. Sd^S} 7 ^ 
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a x. 


-15 

2W0 1°2 7( ' 1 ~ (> 2 


» 8 f 


Jx 


" 5 “ 


103 

2*0^2 (1-P 2 ) 9 2 


a 8 r 
_ 

a x-^ axp 


-io? P 


2 „ 1 «. 2 2 (l-* a ) ^ 


a 8 f 

>X 1 6 SX 2 ? 


1S+90P' 




» 8 r 




- 45 P- 60 p ? 


2 *°; L 6 o 2 4 ( 1 - p2 ) 9 / 2 


3 8 f 


a 7 K *f 
ax 1 «xp 


72 p 2 + 0 + 24 p 4 


S 3, 2x3/2 

2 *o 5 o^(l-p ) / 

X r 


a«r 


X 3 
>X 1 8x 2 


-4*3 p- 6 Op 


2 *o 




a 8 f 

* x ?^ x 2 5 


13 f 9 (V 


2 


x 7 - -2\9/2 

2*o 2 ^ * 


a 8 r 


T 


-103 p 


o r> o 9/ 1 
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3 fl f 


103 


ix 


8 


Q 

2i 0i o 2 9 (1-p^) 


However, since all odd ordered partials of the BVN distribution 
evaluated at the mean are zero, equation (2.4) can be simpli- 
fied as follows 


fCx-pXg) 


2*0^2 (1-p 


*7 


~ jCXt-W 


1 H 1‘ 


2*0, "o 


i 3 ^ 


(l-e 2 ) 


375 


+ ^ ) (x 2 — u ^ ) 


2*° i 2 ° 2 2 ( 1.-0 


2.3/2 ” ^ x 2" u 2 )£ 




+ 25 


2 *°l^°2^ 1 " p2 ^ //2 


- g (x 1 -y 1 )^(x 2 -y 2 ) 


Jo. 


il p p 5/2 

2*o 1 v ? 2 (i-p 2 r / 


+ - 


(2.5) 


From equation (2.5) we observe that ||f(x-pX 2 ) - h(x-pX 2 )ll 
becomes large as (x^,x 2 ) deviates from (y^,^), where h(x^,Xp) 
are the first 25 terms in (2.5) and || • || is some distance 
function. For this reason 

J f(x) - h(x) g(x) dx 
A 

may not be bounded, especially for large region A. However, 
if the regular part h(x) is not removed, the convergence would 
be very slow, l’o accelerate the convergence and allow for 
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integration over large region:;, the PVN program divides the 
original integration region into four rectangular regions 
and integrates each region separately. The program divides 
the four regions as follows. 


Let <_ Xj £ u-^ and Lp <, Xp <Up be the integration 
region. When divided into the four desired regions 
this becomes 


C l 1, u 2 ) 


(l 1 ,l 2 ) 





(Ui,L 2 ) 


Ik 1 u. 

Region 1 limits are <_ <_ — ^5— 

LfiUi 

Region 2 limits are — p — 1 . x i — u 1*’ 

Li +u« 

Region 3 limits are ^ ^ — 

Li+^i 

Region 4 limits are — p — < 


i H U/, 

; L 0 <_ Xp <_ ■ 

i t.i ^ 

^2 - x 2 - 2 
Lp+Up 

; —2— 1 x 2 < u 2 

Lp-fUp 

; 1 x 2 1 Up 


After obtaining the approximate integral for each region the 
results are then added together for the final answer. The final 
standard error is computed as the average of the standard errors 
corresponding to the four regions. 

Since it is difficult to determine if var(h) < 2 cov(f,h), 
the BVN program is currently set up to integrate both the BVN 
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function and the BVN function after extraction of the regular 
part. Convergence is currently checked by computing the 
estimated standard error of 0 after every 1000 random sampler,. 

There are six input items. Those are the means, 
(w 1 ,y 2 ), the standard deviations, the correlation p , 

the maximum standard error, starting value for random number 
generation (odd integer 15 ) , and the limits of integration. 

The estimates for each of the four regions are outputed along 
with their estimated standard error. If the regular part is 
removed, the correlation between f-h and his output. The 
output also indicates whether or not the regular part has been 
removed. Finally, the sum of the values obtained by integrat- 
ing over each of the four regions is displayed as the final 
answer. 


2.3 Specific Examples 

This section presents the output of four examples along 
with the correct answers Pearson (1931). The four integrals 
chosen are 


1 * / 0 / 0 I £» l)dx 1 dx 2 


where E = 


1 .5 

.5 1 


2. / / f(x | o, E )dx n dx. 

$ 1 ~ “ 1 ‘ 


where E = 


1 -.5 

-.5 1 


y 
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5. / / f(x I o, 2 )dx 1 dxp 

o o “ ” 


where £ » 


1 -.75 

-.75 1 


4. 


/ / I £i 
tf l 


£ )dxjdx2 


where Z = 


1 .75 

.75 1 


The results of the BVN program are given in the Tables 
(1-4). 




30 


THE RESPECTIVE MEANS ARE 0.0 0.0 

THE RESPECTIVE STANDARD DEVIATIONS ARE 1.00000000 1.00000000 

THE CORRELATION IS 0.50000000 

THE MAXIMUM ERROR ALLOWED IS 0.00300000 

THE UPPER BOUNDS ARE 4.00000000 4.00000000 

THE LOWER BOUNDS ARE 0.0 0.0 

AN APPROXIMATION FOR THE 1 REGION 

THE VALUE IS 0.2930342318 WITH A STANDARD ERROR OF 0.0025139714 
AND A CORRELATION OF 0.6115916511 

THE REGULAR PART IS POSITIVELY CORRELATED WITH THE INTEGRAL AND THUS 
EXTRACTED 

AN APPROXIMATION FOR THE 2 REGION 

THE VALUE IS 0.0161355461 WITH STANDARD ERROR OF 0.0006924657 
THE REGULAR PART IS NOT REMOVED 

AN APPROXIMATION FOR THE 3 REGION 

THE VALUE IS 0.0164091896 WITH STANDARD ERROR OF 0.0007069048 
THE REGULAR PART IS NOT REMOVED 

AN APPROXIMATION FOR THE 4 REGION 

THE VALUE IS 0.0040517887 WITH STANDARD ERROR OF 0.0002146261 
THE REGUIAR PART IS NOT REMOVED 

THE TOTAL PROBABILITY IS 0.32963076 

WITH A STANDARD ERROR OF 0.00103199 

The correct answer is .33333 


TABLE 1 



THE RESPECTIVE MEANS ARE 0.0 0.0 

THE RESPECTIVE STANDARD DEVIATIONS ARE 1.00000000 1.00000000 

THE CORRELATION IS -0.50000000 

THE MAXIMUM ERROR ALLOWED IS 0.00300000 

THE UPPER BOUNDS ARE 4.00000000 4.00000000 

THE LOWER BOUNDS ARE 0.50000000 1.00000000 

AN APPROXIMATION FOR THE 1 REGION 

THE VALUE IS 0.0111994202 WITH STANDARD ERROR OF 0.0006024142 
THE REGULAR PART IS NOT REMOVED 

AN APPROXIMATION FOR THE 2 REGION 

THE VALUE IS 0.0000608119 WITH STANDARD ERROR OF 0.0000054074 
THE REGULAR PART IS NOT REMOVED 

AN APPROXIMATION FOR THE 5 REGION 

THE VALUE IS 0.0000904058 WITH STANDARD ERROR OF 0.0000072492 
THE REGULAR PART IS NOT REMOVED 

AN APPROXIMATION FOR THE 4 REGION 

THE VALUE IS 0.0000000995 WITH STANDARD ERROR OF 0.0000000122 
THE REGULAR PART IS NOT REMOVED 

THE TOTAL PROBABILITY IS 0.01135074 

WITH A STANDARD ERROR OF 0.00015377 

The correct answer is .012447 


TABLE 2 
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THE RESPECTIVE MEANS ARE 0.0 0.0 

THE RESPECTIVE STANDARD DEVIATIONS ARE 1.00000000 1.00000000 

THE CORRELATION IS -0.75000000 

THE MAXIMUM ERROR ALLOWED IS 0.00300000 

THE UPPER BOUNDS ARE 4.00000000 4.00000000 

THE LOWER BOUNDS ARE 0.0 0.0 

AN APPROXIMATION FOR THE 1 REGION 

THE VALUE IS 0.1118712551 WITH STANDARD ERROR OP 0.0028780424 
THE REGULAR PART IS NOT REMOVED 

AN APPROXIMATION FOR THE 2 REGION 

THE VALUE IS 0.0001379567 WITH STANDARD ERROR OF 0.0000210493 
THE REGULAR PART IS NOT REMOVED 

AN APPROXIMATION FOR THE 3 REGION 

THE VALUE IS 0. 0001607447 WITH STANDARD ERROR OF 0.0000219862 

THE REGULAR PART IS NOT REMOVED 

AN APPROXIMATION FOR THE 4 REGION 

THE VALUE IS 0.0000000005 WITH STANDARD ERROR OF 0.0000000001 
THE REGULAR PART IS NOT REMOVED 

THE TOTAL PROBABILITY IS 0.11216996 

WITH A STANDARD ERROR OF 0.00073027 

The correct answer is .115027 

TABLE 3 . 
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THE RESPECTIVE MEANS ARE 0.0 0.0 

THE RESPECTIVE STANDARD DEVIATIONS ARE 1.00000000 1.00000000 

THE CORRELATION IS 0.75000000 

THE MAXIMUM ERROR ALLOWED IS 0.00500000 

THE UPPER BOUNDS ARE 4-. 00000000 4.00000000 

THE LOWER BOUNDS ARE 0.50000000 1.00000000 

AN APPROXIMATION FOR THE 1 REGION 

THE VALUE IS 0.115327438? WITH STANDARD ERROR OF 0.0027573635 
THE REGULAR PART IS NOT REMOVED 

AN APPROXIMATION FOR THE 2 REGION 

THE VALUE IS 0.0084165793 WITH STANDARD ERROR OF 0.0003595563 
THE REGULAR PART IS NOT REMOVED 

AN APPROXIMATION FOR THE 3 REGION 

THE VALUE IS 0.0033200334 WITH STANDARD ERROR OF 0.0001715015 
THE REGULAR PART IS NOT REMOVED 

AN APPROXIMATION FOR THE 4 REGION 

THE VALUE IS 0.0027903045 WITH STANDARD ERROR OF 0.0003249913 
THE REGULAR PART IS NOT REMOVED 

THE TOTAL PROBABILITY IS 0.12785436 

WITH A STANDARD ERROR OF 0.00085585 

The correct answer is .128133 


TABLE 4 


III. APPROXIMATION 


The introduction briefly presented the reason why 
the Dirchlet distribution was not applicable in the multi- 
variate case. As the beta distribution seemed firmly 
established as a proper model in the univariate case, it 
seemed more reasonable to build a prediction process utilis- 
ing the beta distribution than to seek a new model applicable 
to both univariate and multivariate cases. This led to the 
BVN distribution. 

The reason why the Dirchlet would not work was the 
theoretical requirement of a negative covariance between the 
variables — a situation not frequently encountered in most 
applications. However, the BVN distribution imposes fewer con- 
straints on the value of the covariance. Also, the normal dist- 
ribution has been shown to yield excellent approximations for 
"tail" probabilities in the univariate beta case (See Peizer and 
Pratt, 1968, pg. 1418). Also, the normal approximation exists 
for the beta probabilities over any interval. If the covariance 
(or correlation) is thought of as effecting an increase or de- 
crease in probabilities (compared with uncorrelated probabilities) 
rather than depicting the underlying association between the 
variables, then one should be able to determine this effect 
using either the approximations to the beta probabilities or 
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the beta probabilities themselves. The only reason why a 
bivariate model is required is because we know cloud cover 
frequencies at the sites are related. Otherwise an assumption 
of independence would allow one to compute the joint probabili- 
ties via a direct multiplication of the univariate beta proba- 
bilities. 

Finally, it is important to stress that the BVN, as 
we us, it is only a mechanism to calculate probabilities. 

In conversations with MSFC personnel it was noted that some 
persons in the meteorological profession had proposed the 
normal distribution as a model to describe cloud cover 
frequencies. Such a model may or may not be plausible and. 
we did not investigate it. The beta model serves as the 
basis for our analysis, i.e., we assume the beta model fits 
the data — all we must do is calculate the parameters. Falls 
(1975) did encounter months, time intervals and sites where 
the beta model was not a good fit. It would be proper to 
preface all our remarks and, indeed, the whole report with 
the condition that the beta distribution must yield a good fit 
on the data at hand. However, it is also proper to assert, 
based on proper evidence, that the beta model is always 
adequate, at least for the purposes envisioned. The result 
is the same — situations where the results obtained from 
applying the model differ substantially from empirical results. 
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3.1 Normal Approximation to the Beta Distribution 

Feizor and Pratt (196B) chow that the tail probabili- 
ties for a wide range of distributions can be approximated 
using a normal distribution. Much of the article is not 
germane to our discussion and will not be discussed. However, 
it is informative to trace their procedure for approximating 
the univariate beta distribution. 

The density function for the beta distribution is 

given by 


h(x-,a,e) = 


r (a f b) 
r(a) r(s) 


x a “ 1 (l-x) 8 " 1 ; 0<x<l,a, 


8 > 0. (3.1) 


To approximate the probability that 0<x<x , i.e. 


Pr { x^x Q ) = / 


h(x:o , 8)dx 


calculate the quantities 

d 1 = ( a +8 - 2/3) x Q - (a- 1/3) 


and 


z ** 


1-x^ x -.9 

d 2 , dl , .02 ( / - + , 


5 , l 2( a -t-8-l) — c (e __ 

| 8-.5-(o +8-1) (l-x o ) | 6 (a+ 8 _^j _2 


5)Lof! C, V-iju^To] + 


(a - . 5 ) Log SL1L-? 

[a+ 8 - 1 ] x 


1/2 


(3.2) 


The approximate probability is given by 


,2 


p = r — e-'- y dy. 
^ 2 » 
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Of course, should you desire l;o have a right tail probability. 
The approximate value for the right tail probability 
is 


* - C jk 

The error in these approximations is less than .01 if a, 8 > 1 
and less than .001 if o , 8 > 2. It also follows that 
(x Q < x < x^} can be approximated as 


jC a X 

P { y. <x<x, } - / h(x;° , 8 )dx 1 -/ °h(x:«,B)dx -f h(x;<*, 8 )dx 

TO 1 A V 


or 


! _j z o _i_ _( z i _1_ „-*y 2 


1. '?*■ 


e ' w dy -/ e” Ajr dy =/ 

z, /2* 


o 


dy. 


However, the error is potentially doubled for this case. 

The approximation is not valid for a, 8 < .5 which 
implies the data must be highly U-shapcd for the approximation 
to fail. This could further restrict the applicability to 
some locations and for some seasons. However, Falls has shown 
that this situation is infrequent. 


3.2 The Bivariate Case 

Assuming that x and y are beta distributed, 
x Q < x < Xp y 0 1 y i can be approximated by 


*i . y i 


/. /. 


f( V*y )a V'y • 


(5 


where f (z ,z ) is the BVN distribution defined in equation (2.2) 
x y 


38 


IV. THE APPROXIMATION PROGRAM 

In order to use the BVN approximation, a computer 
program was developed to convert raw data and desired beta 
intervals into the z-values and correlations (BVN program 
inputs). This program takes raw data and calculates means, 
variances, correlations and estimated beta parameters for both 
raw and categorized data. Then for each inputed beta interval 
value (lower and upper values for each variate) it calculates 
a corresponding z-value. 

Two aspects of the program need explanation. Tne 
formulas in Section 3.1 are not defined for the beta values of 
0 or 1. Consequently, the program cannot handle such values. 

For this reason, 0 or 1 must be inputed as 0+e or 1 - r where e > o 
is some arbitrary real number. Likewise -4 is used for - 00 , +4 
for + 00 in the BVN program. 

Since the approximation fails if <* » 6 _< the 

program resets the parameters to .51 and prints a notice to 
the user if the estimated beta parameter value falls below .5* 

It is then left to the user to decide whether or not he wants 
to use this acknowledged poor approximation. 

The beta parameters are estimated using the method of 
moments as described by Hahn and Shapiro (1967, pg. 95 ). The 
estimated beta parameters for the original data are 

7(l-7)-s 2 ] 


B = [ 

A . IS. 

1-7 


39 


_ p 

where X and S are the sample mean and variance. 

A frequency table for both original and category data 
is given in order to compute the empirical probabilities which 
are used to check the corresponding approximate BVN probabilities. 


V. DATA 


The data used in this study was compiled by ESSA, 
National Weather Records Center, Asheville, North Carolina 
and was provided to the authors by Organization ES-42, 
Marshall Space Flight Center, Alabama. The sites selected 
were Fort Worth and Houston, Texas. Daily records (January 
1971 to December 1975) on cloud cover, measured in tenths, 
were recorded every third hour. 

The data was grouped into the categories shown in 
Table 5 (Fall 1975). 


Table 5 

Cloud Cover Categories 
Category Tenths 

1 0 

2 1 , 2,3 

3 4,5 

4 6 , 7 , 8, 9 

5 10 


Since Falls (1971) demonstrated that the beta distribution 
adequately describes variation in categorial data, our primary 
investigation was restricted to categorical data. However, 
the approximation program is not restricted to categorical 
data. 
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VI . EXAMPLES 

A complete set of probabilities (25 values) have 
been calculated for the Port Worth 9 a.m. and Port Worth 
5 p.m. combination. These values are presented in Figure 1. 

Each of the five portions of figure represents a category 
level for 9 a.m. and the abseissas represent the categories 
for 3 p.m. Table 6 presents a portion of the approximation 
program and Table 7 gives the corresponding BVN computations. 

Figure 1 values were determined based on observed and 
expected frequencies for 5 years (155 values). As can be 
noted, the agreement is quite satisfactory with a couple of 
exceptions. Values for Category 1 for 9 a.m. and Category 2 
for 3 p.m. shows a wide divergence. Also the five values 
predicted for 3 p.m. and Category 4 for 9 a.m. show substantial 
disagreement. 

These discrepencies between observed and predicted 
values can be explained by analyzing how well the beta model 
describes univariate cloud cover in the various data sets. 

Prom Table 6 the category frequencies for Site 1 ( 9 a.m.) 
are 46, 29, 20, 39, 21 respectively and the estimated beta 
parameters are . 862646 and 1.06241. These parameters are for 
a very U-shaped density which decreases as x 1. Consequently, 
the fitted distribution does not reflect the variation in these 


Frequency - 9am. Frequency - 9a.m. Frequency - 9a.m 


2 3 4 

Categories - 3 p.m. 

CATEGORY 3 


2 3 4 

Categories - 3 p.m. 

CATEGORY 4 



2 3 4 

Catogorioo - 3 p.m. 
CATEGORY 5 


9a.m. Cloud Cover Mean s .443387 
Standard Deviation -- .P.90846 
3p.m. Cloud Cover Mean = .510323 
Standard ue'/iation= .24098/ 
Correlation = .570933 


Beta Parameters for 9a.m. are .862646 1.061241 

Beta Parameters for 3p.m. are 1.685583 1.617392 

FIGURE 1. 

OBSERVED AND PREDICTED FREQUENCIES FOR FORT WORTH AT 9 A.H. 
AND FORT WORTH AT 5 P.M. BASED ON JULY DATA FOR 1971-75- 
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data for categories 1 and 4, which is reflected in the approx- 
imate probability. 

Some additional comments are necessary. First it is 
important to note that we have only 155 data points and more 
data would, in most such cases, give better fit to the true 
distribution hence a better approximation. Secondly, this 
problem is not restricted to this one isolated case. Based 
upon our analyses, we feel that the substantial disagreement 
between observed and predicted probabilities were based upon 
the inadequacy of the beta distribution. It does not seem 
likely that large errors will occur because of this condition 
but if the parameter values are low the approximation error 
could contribute substantially to the disagreement between the 
values. Thirdly, it must be noted that Figure 1 is based upon 
integration limits (determined by the transformation from 
categories to the (0,1) interval) that should give the best 
results. The category values 1, 2, 5, 4, 5 are transformed 
to .1, .3, .5, .7, .9 respectively. The corresponding limits 
of integration are found in Table 6. 

Table 6 

Integration Limits 

Category Integration Limits Midpoint 


1 .01 to .2 .1 

2 .2 to .4 .3 

5 .4 to .6 .5 

4 .6 to .8 .7 


5 .8 to .99 .9 
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The values in Table 6 are the usual "continuity" corrections 
for approximating probabilities for discrete variables. It 
must be noted that the intervals selected will not always 
reflect the underlying situation and hence could contribute 
to the differences in values. However, if the above limits 
are a source of error then its effect will be minor compared 
with the other errors and its effect will decrease over wider 
intervals. 

As noted, we have elected to use categorical data 
throughout the analyses. However, one might consider using 
the original data in that the beta model might actually fit 
whereas the categorical fit was inadequate. Another reason 
for using the original data is the greater flexibility in 
selecting the integration limits which can be made to closely 
agree with the original situation (cloud cover measured in 
tenths). 

6.2 Application of the Programs 

The approximation programs must be run to obtain the 
approximate integration limits used in integrating the BVN 
distribution. The input needed for this program consists of 
two parts. The first part consists of the raw data (read 
pairwise with the first value corresponding to the first site 
and the second value corresponding to the second site or the 
data can represent one site at two different times). The 
second part consists of the input ed boundary numbers for the 


regions to be integrated. Before continuing one should 
inspect the outputed beta parameters and corresponding 
frequency tables. If the estimated beta parameters are 
significantly less than .5, then one must proceed with 
caution since the calculated integration limits are probably 
unreliable (for reason explained previously). 

The outputed correlations and the integration limits 
are then used as inputs into the BVN program. Note that 
since the approximated integration limits pertain only to 
the standard normal distribution, the mean vector will be 
(0,0) and the standard deviation will be (1,1). The main 
output of the BVN program is the total probability. This 
value represents the approximate probability of a specified 
category or categories at Site 1 intersected with a specified 
category or categories at site 2. 

For example, Table 7 lists the output of the approxi- 
mation program for the percent of cloud cover over Fort Worth, 
Texas, at 9 a.m. and 3 p.m. during the month of July (1971- 
1975)* Since the beta parameters for the original data is 
significantly less than .5, we decided to work with the 
category data. The category data z's are the approximate 
integration limits corresponding to category 1 at 9 a.m. and 
category 1 at 3 p.m. These values were then used as input 
for the BVN program along with the correlation of .57. The 
output of the BVN program is found in Table 8. The total 
probability of having cloud cover in category 1, (i.e. 
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essentially no cloud cover) at 9 a. in. and of having cloud 
cover in category 1 at 3 p.m. during July at Port Worth is 
shown to be approximately .063. Whereas the empirical value, 
found in the category frequency table, is = .0645. 
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JULY 71-75 r T . 4fJMT H 9 A.M. AND 3 P.M. 
♦***♦ RfcSULTS USING ORIGINAL DATA ***** 


FREQUENCY TABLE » 155 VALUES 



10 

3 

9 

11 

9 

1 

0 

3 

3 

0 

0 


1 

0 

1 

2 

4 

2 

0 

2 

0 

0 

1 


0 

0 

0 

1 

2 

2 

0 

2 

1 

0 

1 


0 

0 

1 

2 

0 

1 

0 

3 

1 

2 

2 

0 


2 

2 

0 

I 

4 

2 

1 

0 

0 

1 


1 

0 

0 

1 

0 

1 

0 

1 

<j 

0 

0 


0 

0 

0 

2 

0 

1 

0 

2 

0 

0 

1 


0 

Q 

2 

5 

1 

2 

0 

0 

2 

1 

0 


0 

0 

0 

l 

0 

5 

0 

2 

i 

0 

4 


0 

0 

c 

0 

0 

0 

2 

1 

1 

3 

5 


0 

0 

0 

0 

1 

3 

0 

0 

** 

5 

8 

MEANS* 



0.4151 


0.5065 







ST. DEV* 



0.3798 


0 . 316 ? 







CORR • 



0 . 64 1 1 









ESTIMATE) 

BETA 

PARAMETERS 









SITE I 



0.2849 


0.3998 







SITE II 



0. 7603 


0.7406 








*♦♦4* RtSULTS USING 

CATEGORIC 

AL DATA *♦•** 

FREUUEN 

CY 

table; 

155 VALUES 

10 

23 

10 

3 

0 

1 

7 

11 

8 

2 

3 

5 

7 

4 

1 

0 

10 

4 

15 

10 

0 

0 

4 

9 

8 

MEANS* 


0. 4494 


0.5103 

ST. DEV- 


0.2908 


0.2410 

CORK « 


0.5739 



SITE I 


0. 8626 


1.3612 

SITE 11 


1.6855 


1.5174 


***** NOTE ***** 


IF A PARAMETER OR PARAMETERS IS LESS THAN OR EQUAL TO .5* THE Z-VALUE IS 
UNDEFINED. FUR FURTHER COMPUTATION THE PARAMETER IS Rf " “ 


tESET TO .51 


FIRST SITE 

UPPER LOWER UPPER 


StC0N0 S1 l5viS* 


INTEGRAL LIMITS 
CATEGORY DATA 2‘$ 
ORIGINAL OATA Z'S 


0. 90030000 
1.02438 124 
0.69792 5 24 


0.20000000 
■0. 64939553 
■0.69792529 


3. 80000003 
1.08935735 
0. 714J007J 


0.20003033 
-I. 17765974 
-0. 75874548 


TABLE 7 
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INPUT PARAMETERS 



MEANS- 

0.3 

0.0 

ST. DEV- 

1.3000 

1.0000 

CORRELATION- 

0.5709 


MAX ERROR- 

0.3030 


UPPER 0QJNOS- 

1.3241 

1.C894 

LOWER 8U JNDS * 

-0.6494 

-1.1777 


***** AP»ROX. FOR REGION NO. - 1 


♦♦♦the regular part has been re moved 

THE VALUE IS 0.1 3776 
ST. ERROR « 0.00001 
CJRR • 0.56849 


♦ *♦** AP»ROX. FOR REGION NO. • 2 **♦♦* 


♦♦♦THE REGULAR 
THE VALUE IS 
ST. ERROR » 
CORR - 


PART HAS 
0.07902 
0.00065 
0. 55691 


BEEN 


REMOVED 


***** AP>ROX. FOR REGION NO. * 


***** 


♦ ♦♦THE REGULAR PART HAS BEEN REMOVED 

THE VALJ= IS U. 12317 

ST. ERRU* » 0.00006 

CORR ■ 0.66100 


***** APPROX. FOR REGION NO. ■ 


***** 


♦♦♦THE REGULAR PART HAS BEEN REMOVED 

THE VALUE IS 0.13429 

ST. ERRuR - 0.00000 

CORR • 0.74074 


THE TOTAL PROBABILITY IS 0.47424 WITH A STANDARD ERROR JF 


0.00018 


TABLE 8 
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APPENDICES 


Appendix A gives a description of the card inputs 
followed by a listing of the BVN program. Appendix B gives 
a similar listing for the approximation program. Both programs 
are written in Fortran and, the approximation program can 
compute 100 individual integration limits in less than a 
minute (on IBM 370 / 155 ). The BVN program also takes less 
than a minute to calculate one total probability. 
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APPENDIX A 

BVN Program - Card Input 


Card 

1 

1-14 



15-28 

Card 

2 

1-14 

15-28 

29-42 

Card 

3 

1-14 

21-25 

Card 

4 

1-14 



15-28 

Card 

5 

1-14 



15-28 


mean of the first variable of the BVN 
distribution 

mean of the second variable 

standard deviation of the first variable 
standard deviation of the second variable 
correlation between the two variables 


maximum standard error allowed 
odd, five digit random integer 

upper integration limit for the first 
variable 

lower integration limit for the first 
variable 


upper integration limit for the second 
variable 

lower integration limit for the second 
variable 


i 

■1 

i 


| 


y 
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g #♦•#*#♦♦♦♦**♦ ********************C ARD 
8 


input ♦+♦++*♦********* **♦*****♦**♦**♦**♦♦# 


c 

c 

c 

c 

8 

C 

c 

c 

c 

c 

c 


CARO COLS. VARIABLE 


1 

1 

z 

2 
2 
3 

3 

4 
A 

5 
5 


1-14 
15-28 
1-14 
15-28 
29-42 
1-14 
21-25 
1-14 
1 5-28 
1-14 
15-28 

IMPLICIT 

COMMON 


■VAN AT FUST SITE, REAL NUMBER 
■45AM AT SECOND SITE* RtAL NUMBER 
STO.UEV. AT FIRST SITE* REAL NUMBER 
STD. OEV. AT SECOND SITE* REAL NUMBER 
CORRELATION# REAL NUMBtR 
ERROR BOUND 
INTEGER RANDOM NUMBER 
U®PER BOUNO AT FIRST SITE 
LlWF.R BOUND AT FIRST SITE 
UPPER BOUNO AT SECOND SITE 
LOWER BOUND AT SECOND SITE 

REAL*3t A-H*U-Z) 

B, Ml# M2#K0* S IG1, SIG2 


R t A l. * 8 C(Hl/a*0./»X(4)#Nr#Ml#M2»8<8»R)#QC4)/4*0./#LUU4),UPl(4)»LJ 
C 2 ( 4 ) » UP 2 <4 ) * LI » L2 
INTEGER RAND 

15 READ! 5* 1 *6 *10 = 100) M l » M2» S I G 1 # S I G2 # R J# E RR DR # R AND# U 1 » L 1 » U2 . 1 2 

1 FORMAT! 2 FI 4, fl/3F14.«/F 14.8# t>X» I 5/ 2F 1 4. 8/ 2F 1 4 . 8 ) 

WR I TE ( b #?0D ) 

,-OD FORMAT ( * INPUT PAR AMETFRS • # // ) 

WRITtO.2) M1*M2#SIG1* SIG2»R3 

2 FORMAT! 1 MEANS* *# 20X# 2Fl J.4#/» • ST. DE V « * * 1 3 X* 2 F 1 0. 4 » / * 

♦ • CORRELATIONS*# 14X»F10.4I 
WRITE lb# J) ERROR 

3 FORMAT!* MAX E RR DR * * # 1 6X . F 1 0 . 4 ) 

WKI T E ( b » 4 J J1*U2*L1»L2 

4 FORMAT ( ' UP 0 ER BJJNOS= * * 1 3X# 2F 10. 4# /* 

* • LOWER BOUNDS*' *13X»2F10.4» //) 

JP1I1)*< L1+J1I/2D0 

UP 1 ( 2 ) * J l 
JP1I3) = JP) () ) 

UP1 1 4 ) * U 1 

UP2( 1 >»( L2 + U2) /2D0 
JP2(2)*JP2 ( l ) 

UP213 )=U2 
JP2 C 4 ) = U 2 
L01I1I--L 1 

L01(2>*( L1+U1I/2D3 
LOU 3 ) *t 1 
L01I4I-L Cl (?) 

L 02 I 1 )*L 2 
L02 ( 2 ) *L 2 

L 02 ( 3 ) * ( L2+U2J/2D0 
L02 l 4 ) *L C2 ( 3 ) 

B ( 2 * 1 ) * - l.ODO 

B ( 2 » 3 ) *- 1 . ODO 

3(2#2)=R0 

B(4»1)=3D0 

B(4*5)*3DO 

B ( 4 * 2 ) »- 3U C* RO 

B(4* 4)*- 3D 0* RU 

B ( 4# 3 ) =2 DO* RO* *2 <-100 

B(6#l l«-15DD 
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B( b» 7 ) ■- 1503 

3(b» 2)»l 5004*0 

B( b» 6) * l 500*R0 

3lo»3 >«- 3u0- 12D04RU*42 

d(o#4 )*RD0MiH6D04R'J443 

3(6.5)*- 3DC- l?00*RU**2 

d(rt»l> = 105D3 

3(rt.2)»(-10500)4R3 

3 ( S » 3 )*l 5D0490004RQ442 

3(8»4)*-45D34RU-bOD04RO**3 

3( 3. 51*720 C4RU442+24304R0444+9D0 

B(3.b)*B(8#4) 

3(3.7)*3 I d > 3 ) 

3(rt,3)*a(8.2) 

3 (3. 3) =3 (8.1) 

;U|-.| UGD0)/UD04 3. 1415926535900*SlGl4SlG24ROS0**( .533) ) 
DO 8 J*l.4 


J J * 2 * J 

b :ijji<:i n /<jsu4*j 

$U=330 


St* 330 
JG 16 1 * 1# 4 

16 CALL T AY LUR U ( I )> C. JP1 ( I ) » UP 2 ( 1 ) . L3 1 ( l ) » LQ2 ( I ) . l) 
3U 17 1*1# 4 

18 FORMAT?//?' ***** APPROX. FOR REGION NJ. **#14# • 
REGP=0< I ) 


NT= l 003 
CGF 0= 000 
GSUV0J3 
F SUN* 033 


***** •, // ) 


FGS 3=000 
FSQ = 000 
PRO* < UP l 
00 5 II* 
00 7 J=l 
CALL 
RAN 3 


( I ) -L 0 1 ( I ) )MUP2( I 1 -L 02 ( I 
1# 1000 
>2 

RAN DU ( R AND# I Y, YFL ) 

iy 


) ) 


X(J) =YFL 
X( 1 ) = LOl (1 )+Xm*(U<>l( 
X ( 2 ) = L02 (I )*X(2 )4(UP2I 


I )-L01 ( I ) ) 

HOL * ( -1 . 666/ ( 2D04R0SQ ) ) *< ( X( 1 ) -N1 ) 4 *z / S I 3 1 **Z-Z D04R0* (XU) 
C(X(2)-N2)/(SIGI*SIG2)+(X(2)-M2)**2/SIG2442) 


,Ni )* 


F = C ( l ) 40 EX P ( HOL ) 

F*F*PRO 

FS0*FS3* (F 4* 2 ) 

FSUN=FS J M*F 

CALL TAYL0fi(G»C#X(l)»X(2)»L01(I)»L02(I)»2) 

G=G4 PRO 

GSUN *GS J M+G 

FGSQ=FGSQ+ (( F-G)442i 

COFG=CJF G*F4( F-G) 

6 CONTINJE 

FGSUM=FSUM-GSUM 

FGM = FGS J M/NT 

FVAR=(FSQ-(FSUM 442/NT) ) / ( NT-l.ODO) 
VVAR=(FGS0-FGSUM**2/INT)/ (NT-1. 000 ) 
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10 

a 

12 

m 


m 

14 

li 


.000) 

SF»'0SQRT IFVAR/NT) 

Crt-tSUi/NT 

IF ( SFG-ERR3R) 9,9,10 
1 F I S F-t R ROR ) 12.12. 13 
MT»NT«-l 000 

oo to 5 

CONTINUE 
WR l TE ( b» 1141 

FrjRiAT 1 1 ♦♦♦THE REGULAR PART IS ROT REMOVED 1 ) 

WR I T E I b, 14) FM,SF 

SE«$E+SF 

SU-SJ+FR 

50 TO 17 

FG« e GM*RFGP 

C0F-C0F7US0RTC VVAR*FVAR ) 

$E»SE+SFG 

$U«S'J«-FG 

111) 

♦♦♦THE kEGJcAR PART HAS BEE< RF NOVEJ • ) 
ll) FG,SFG,COF 
THE VALUE l $ 1 , F 10. 5, /» * ST 


HR I r E ( t> » 

FDRNATI • 

VRlTSlb, 

FORI AT ( * 

F0R1AM • 

* • CiiRR 

17 CONTINUE 
SE-SE/400 
WRlTElb. 19) SO, SC 
19 FQRRATI • O’ /• OTHE TOTAL 
COR OF * , F 12.5) 

50 TO 15 
100 STO* 

ENO 


THE VALUE 1 S ' * F 10. 5, /, ' ST. 
» ', 3X, F 1 1 • 5, / ) 


ERR3R 

ERROR 


* » F 1 1 • 51 
' * F 1 1 . 5, / » 


PROBABILITY IS', FIG. 5,' WITH A STANDARD ERR 


?HPi , ?l/R N I 6 im?2.A?a^V l ’ u2 ' Ll ' L2 ' ,NC ' 

C0MN3N B* Ml# 12, RO, S IG1, 3 IG2 
RfcA'.*H C (8) . L1,L2, Ml,N2,3t 8.9) 

J*C( l ) 

IF (INC.EU.l) Q* , «1*IU1-L1)*(U2-L2) 

00 12 *»2, 3,2 
NK = <* 1 

V AR 2» l . 0 DO 
DO 12 J=1,NK 
JJ=J-1 

IF ( J.LE .2) GO TO 17 
J3«J-2 
NVAR2-JJ 
DO 15 L = 1> J 3 
15 NVAR2«NV AR2MJJ-L) 

VAR2«NVA R2 

17 IF IJ-iO 18,19,19 

18 NVAR3*(<-JJ) 

KH*X- J J - 1 

00 lb L= 1, KH 
lb NVAR3*NV AR 3* IK- JJ-L ) 

V AR 3 =NV A R3 
30 TO 20 

19 VAR3-1.0D0 

20 VAR«l.OOO/IVAR2*VAR3) 

IF l INC. EQ.2 ) GO TO 14 

0«OMVAR*IC(K)/(SIG1**(K-JJ )*SI32**JJ) ) 

C*(IU2-M2 )**J-(L2”M2)4*J)*( 1.000/1 (N<-JJ)*J) )4((U1-M1)P*(NK-JJ)- 
C(LI-N1)**(NK-JJ))*B(K»J) ) 

GO TO 12 

14 OOMVAR*(C(K)7(SIGL**(K-JJ)*SIG2**JJ)) 

C*( 1 J2-N2 )**JJ 1M (U 1-MI I44INK-J I ) *B( K, J ) ) 

12 CONTINUE 
RETURN 
END 
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Card 


Card 


Last 

Card 


APPENDIX B 

Approximation Program - Card Input 


Cols. 


1 1-4 

5-80 


2 + 1-76 


1-10 

11-20 

21-50 

31-40 


number of data pairs 
19 pairs of data with each element of 
each pair right justified in a two 
column space; no decimal points 


19 pairs of data with each element of 
each pair right justified in a two 
column space; no decimal points. That 
is the data is read with an 19F2.1 
format . There will be as many cards of 
this type as necessary to punch all data. 


lower integration limit 
upper integration limit 
lower integration limit 
upper integration limit 


for the first site 
for the first site 
for the second site 
for the second site 


noononnnnnnonnivr) 
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t**«**««*«»**t**t***t»****t*****CARO I M PUT* ************** *••••*•••**•****•*♦— 
CARO C31S. VARIABLE 

1 1-83 TITLE 

2 1-5 NJMfttR OF DATA FROM EACH SITE 

2 5-80 ALTERNATING DATA; FIRST SITE CLOUD COVER THEN SECOND SITE CLOUD 

COVER ttOTH AN INTEGER BETWEEN 1 AND 10 IN CONSECUTIVE T*U COLUMN 
S® AC FS . 

3* 1-76 CONTINUE DATA INPUT AS ABOVE 

A 1-13 LOWER INTEGRATION LIMIT FOR SITE ONE 

A 11-23 U»PFR INTEGRATION LIMIT FOR SITE ONE 

A 21-03 LOWER INTEGRATION LIMIT FOR SITE TWO 

A 3 1 - A 3 U^PER INTEGRATION LIMIT FOR SITE TWO 

ALL OATA ON CARD A MUST BE LESS THAN l SINCE W£ ARE DEALING WITH THE BETA 
DISTRIBUTION UNO GREATER THAN ZERO) 

IMPLICIT REAL*8 (A-H#M#0-Z) 

REAL *8 X 1155)* Y( 155) »CX( 155), C YU 55) #MX. AAl 13 ) 

INTEGER F(5»5)/2S*0/»Fa(ll»ll)/l2l*3/ 

RE A3 1 5» 76 ) ( AA( 1 ). I«l, 10) 

WR l T E ( 6» 7o) (AAU)>I*1»10) 

76 FQRNATI1 CA8) 

59 READ! 5, 37»END«1001 N . (X d ) > Y II) * I « l # N ) 

37 FUR1AT ( IA » I 3BF 2 • l )) 

30 60 1-1*5 
30 63 J« 1» 5 
63 F ( 1 » J ) *3 

30 71 1- 1* 11 
30 71 J* l> ll 
71 FOI I* J ) - 0 
30 72 1 — 1 * N 
Nl«XI I ) * IODO + 1. IDO 
N2«Y( I)*10D3M. IDO 
7 2 FJHUN2 )>F ]( N1 »N2 ) + 1 
CALL CAT ( X » C X » N ) 

CALL catiy»:y»n) 

JU 5L I = l* N 
N 1 * C X ( l) 

N2 * C Y ( I ) 

51 F ( N 1 » N 2 ) *FtNl»N2)+13Q 
WR I T E ( 6» 200) 

203 FORIAK//*' ***** RESULTS USING ORIGINAL OATA «****•, //) 

WR I Tc ( b* 73) N 

73 FOR NAT <'0'.22X.' FREQUENCY TABLE ; • » I 5* 2 X» ' VALUES • » // ) 

30 73 1- 1* 11 

73 WRITEI6* 52 ) ( F (J ( I » J ) / J = 1 » 1 1 ) 

52 F0R1AT ( 10X. 11 l 5) 

30 55 I » 1* N 

CXUIMCXI D-.5D0) /5D0 
55 C Y ( I) «<CY( I)-. 500) /5D0 

CALL STA T(X»Y»MX»MY*$VX»SVY»R01»N) 

CALL STAT ( C X, C Y, MC X» MC Y, S VC X , S VC Y , R 02 » N ) 

SIG1-3S0RTISVX) 

SIG2=0S3RT(SVY) 

S IGCUDSORTtSVCX) 

SIGC2-DS0RTISVCY) 

WRITE (b;48) MX.MY/S IGUS1G2. RUL 


”• MV...M.ZMO.W. 

1XJ-SVX) 


201 

50 

2 02 
5* 


Hi* If IDO-MX) /SVX) *INX*t IDO- 
Al«MX«"U)/( IDG-MX) 

J2«l! I JJ-MYI/SVY) *(MY*< ID0-NY)-S¥Y) 
A2»!MYM2)/I 100-MY) 

WRITE !b# 201) 

L- S r I M A T E 0 JET A PARAMETERS • 
50) Al#Bl»A2»82 


format ( • 

WRITE !o. 
FORMAT! • 
WRITE! o# 


// I 


l •» 10X.2F10.4# /# • SITE I 1 * #9X# 2F10.4#// ) 


US 1 NO CATEGORICAL jaTA *****•#//) 


53 


203 

75 


55 

55 


SITE 
202 ) 

FORMAT!//#' ***** RESULTS 
FORMAT ( 1 0‘ # 10X# 'FREJUENCY 
WRITE I 6# 54, ) N 
JO 53 I * I# 5 

WR I Tc I b» 52 I (F( I# J ) # J* 1# 5) 

WR I T E 1 1>» 48 ) MCX.MC Y. S 1 GC l . S I GC 2# R 02 
HBl»( ( I0 0-MCX) /SVC X) *1 MC X* I 1 DG-NC X ) - S VC X ) 

MAI « ( MCX *HB l ) / ( 1 DO-MCX ) 

M82 * I (10C-MCY)/SVCY)*!MCY*( 1 DC-MC Y ) - S VC Y ) 
rlA2*( MCY*Hi? ) /( 1U0-MCY) 

WRITEI6.50) HA 1 » HB 1 » HA2 » HB2 
WRITE! 6. 203) 

FORMAT!//#' ***** NJTE **♦**•) 

WRI TE ! 6 # 75 I 

FORMAT!' ' » ' l F A PARAMETER OR PARAMETERS IS LESS THAN HR EJUAL TO 
C .5# THE Z-/ALJE IS'/' UNDEFINED. FOR FURTHER COMPUTATION TriE PARA 
CMETEI IS RESET Til .51*1 
2 E AO ! 5»56#tN0=10J) Y1#Y3#Y2#Y4 
c 0 R M A T ( 8E1 0.01 
ALL CALZ(H3 1#HAl#Yl#N#ZCl) 

CAL Z( HJ2#HA2» Y2#N# ZC2I 
CAl U HJ i#HAl# Y3#N# ZL3) 

CAL Z( H32#HA2»Y4#N#ZL4) 

CAL 7 1 8 l # A 1 » Y 1 # N » Z 1 I 
CAL ZIB2#A2#Y2#N»Z2) 

CALZ (5i»Al#Y3»N#ZLl) 

CALZ!82#A?#Y4.»N#ZL2) 


call 

CALL 
CALL 
CAL L 
CALL 
CALL 
CALL 
WRITEIb# 74) 


74 FORMAT! ' O'# 3 l X# ' F IRST S I TE • # 20X. ' SEC JUO S I TE » / 26X. ' UP PE R ' # 1 OX# • L UW 


70 FORMAT!* 

CRY JATA 
C2« F14.S# 

WRI TE (b» 70 ) 
68 SO TO 58 
100 STOP 
END 


8, IX# E 14. 8# IX >/lX# 'CAT EGO 


C E R ' » 10X# ' U P a E R ' »10X# 'LOWER' ) 

O' #2X» • INTEGRAL LIMITS 1 » 1X»2( F14, 

Z' 'S' #1X»2IF14.8#1X#E14.8#IX)/1X# 'ORIGINAL DATA Z''S'#IX# 
1X»=I4.8# IX) ) 

Y3»Y1#Y4#Y2#ZL3»ZC1»ZL4»ZC2»ZL1#Z1#ZL2#Z2 
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*7 


SUB33UT1NE S TAT ( X / Y , M l , M2/ SV ARX » $ VARY# RO. N ) 
IMPLICIT REAL*8 4A-H/M/U-Z) 

X(155l»Yllj5l 


REALMS 
SxODO 
SYOOO 
>XS *000 
SYS-3JU 
SXYOOD 
30 *7 I « 1# N 
SX»SX*X( 1 ) 

SY»SY*Y( I ) 

SXS-SX5-*- I X< I 11**2 
SYS-SYS*- CY« I ) )**2 
SXY«SXY*X( I I * Y( I ) 

H»SX/ <4 

> 1 ?» s y/n 

5>VA<X-(SX5-'(*M1**2 )/ (N-1D0) 
5VARY»(SYS-( $f**2/NI )/ ( N- 100 ) 
R0«( GXY-N*;U*M?I/(N-130| 
RO-RO/OSORU $VARX*$VARY) 

RET J<N 
END 


5y3igdTI NE CAlZtBl.Al.Y.M.Z) 

IMPLICIT RcU*8 (A-H/M/O-Z) 

IF (Al.LE.53-1) A l * 5 1 0-2 

IF (tJl.Lt.53-l) 81«513-2 

SX1 « B1-. 500 

SX2 * Al - . 5D0 

SXN-A1 + B 1-1-30 

P* L 33- Y 

/.?!T ( §?'. ,| t* 33333;,3;i00,<lY - ul '* 3 333333 300)*20-2*(Y/Bl-P/AU(Y- 

C I A l *3 I ) 1 

DA* DABS ( $X 1 - $XN *P ) 

3L5* JL J3 (SXl /( SXN*P) ) 

3LT* JLJG (SX2/4 3XN*Y) ) 

Z<O2/0A*0S0Rr(l2CO*SXN/(63O*SXN«-100)*( SX l *0L S + S X2*0LT ) I 

R E T J R N 

ENO 


*v : ,, 


mimtu 

30 38 I* 1/ N 
IF (XU)-iDO) 39/ A 0» AO 
AO H( I ) *530 
GO TO 38 

39 IF (XU )-.630)Al/A2» A2 
A2 HMI-AOO 
G G TO 3 5 

Al IF (XU I-.A30) A3» AA» AA 
AA H( 11*300 
GO T D 3 8 

A3 t F (XU )-. 130) A5/ A6/A6 
Ab. H( I) = 230 
GO TO 38 
A 5 H( I » * l DO 
38 CONTINUE 
RETURN 
END 


y 


. 500) / 
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A PROGRAM TO COMPUTE CONDITIONAL BIVARIATE 


NORMAL PARAMETERS 


Summary 


This report derives the conditional bivariate normal 
parameters from an original quadravariate distribution. The 
paper presents the theory and appended is a computor program 
developed to give numerical results. An example is presented 
in the paper. 


I. INTRODUCTION 

This report presents a sketch of the theory and a computer 
program designed to calculate the bivariate normal conditional 
distribution derived from the quadravariate normal distribution. 
The required computor inputs arc described and an e v ..mplo is 
presented. The computer program is appended. 

Theory 


The general multivariate normal distribution has the density 

' r-1 


f (t j > * • • • 


1 


(2*) k/2 m 


I ' l TT f/ T exp I" Ql- m ) } (D 


where p_ « (p , p 2 p^) , the vector of mean values and 


! o n o J2 ...o lk 


°21 °22 •'•° 2 k 


°kl °k2 ’°kk 
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A property of the multivariate normal distribution is that 
marginal and conditional distributions are also normally distributed. 
The general expression for these distributions are found often in the 
literature Csec Morrison (1967)] . We shall confine remarks here to 
the specific case. 
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then letting 
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Where 


and 
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2»! 

r*l 5 

r 
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we have 


* f * _ 1 

X , X * 


a. • h * E 12 E 22 ‘*2 


( 2 ) 


(*) 


(4) 


Computation of the parameters for this conditional distribution really 

* * 

reduces to computation of the quantities £ and . Carefully note 

* r -i' 

that the value of u includes values of x_ = Lx_, x.J that must be 
— —2 3 4 

* 

specified before numerical values for _v can be calculated. 

+ 

Even for this rather easy case the actual expressions for £ 

* 

and £ and therefore for the quadratic form in (2) are very romplicatcd 
algebraically. They arc, however, very amenable to numerical compu- 
tation via computer. The least complicated for the expressions is 

* 

that for jj and the actual form is given below (letting 0 34 bo 43 f° T 
convenience ) • 
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V U, ’l*°4 • ,, 14 , \W Mx r M 3 ) + (< ’l4^3-°13 < i4 )tX 4-‘ , 4 )} ' 

V* ,0 2S 0 44'°24°33 ) + *°24 ^5“°25°.W * 4* u 1 ) * ,(o 33 0 44“°34 ) 


_i * 

The matrix triple product ^ 22 ^21 wa * tes * a complicated 
expression and this, of course, causes ( E ) and, therefore, the 
quadratic form in (2) to be almost incomprehensible in an expanded form. 


Computer Program and Required Inputs 


The computer program is written to accept quadravariate data 

and return the conditional bivariate parameter. The conditional variance- 

covariance matrix and the associated standard deviations and correlations 

arc initially calculated and printed. The program is designed to take 

as many pairs of "conditioning values" of x^ and x^ as desired and print 

* 

out both the values of x^ and x^ plus the associated values of u_ . 

Example: The following data was input to the program 


p * 

[21.58, 

-.04, 

, 43.35, 

1.25]' 

*°i7 * 

11.03, 

M 

c* 

Q. 

.0503, 

P n - .7382, 

So 2 2 “ 

11.52, 

P 23“ 

1614, 

P 24 - .8134 

• / °I7 ■ 

15.47, 

P 34‘ 

.1524 



14.59 




[Xj.I 

< 4 3 - C4: 

*.35, 

0] 



Attached as Appendix I is the output giving the calculated parameter*' 
for the bivariate conditional. Note carefully that the standard 
deviations and correlations are printed in matrix form for convenience- 
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not to he confused with the variance-covAriancc matrix printed above it. 
Below the standard deviation and correlation matrix the values condi- 
tioned on and the resulting conditional means are printed. The original 
inputs and matrices will be printed only once hut the values conditioned 
on, followed by the conditional means calculated using those values, 
will he repeated for each set of conditioning values read in. 

Input to the program consists of the followin'; cards: 

Card 1 The 4 means for the quadravariate normal in 4F10.4 
Format . 

Card 2 Standard deviation for variable 1 followed by 

correlations for variables 162, If, 3, and 163 in 
4F10.4 Format 

Card 3 Standard deviation for variable 2 followed by 

correlations for variables 263 and 264 in 3F10.4 
Format . 

Card 4 Standard deviation for variable 3 followed by 

correlation between variables 364 in 2F10.4 Format. 

Card 5 Standard deviation for variable 4 in F10.4 format. 

Card 6 Number of sets of x^. x. values to be conditioned 
on in 12 Format. 

Card 7 l St set of x_, x, values to be conditioned on 

3 4 

Card 8 2 nd set of " " " " " " " 

•• i, «• it I, », », it ii «i 

" etc. 

The source deck listing is given in Appendix II. 

References 

Morrison, D. F. (1967). Multivariate Statistical Methods . Wiley, 

N. York. 


«* * 


62 


APPENDIX I 


MEANS VECTOR • 21.5800 -0.4000 43.3300 1.2300 

V AR I Af 1C E - CO V AR I ANC E MATRIX 


121.6009 

G.3914 

125.9621 

-3.2025 

CONI). VAR. COV. MATRIX 

53.17073 

4.C3C23 

SU&CORR. MATRIX 

7.29244 

0.00022 

values conditioned nr; 

CONDITIONAL MEANS 


6.3014 125.0621 
132.7104 20.7630 
28.7630 230.3200 
136.7336 34.3*70 


4.03024 

44.71625 


0.00022 

C.6C7C2 

43.3500 -0.0070 

21.7001 -C.CC70 


-3.2025 

136.7336 

34.3070 

212.0631 
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APPENDIX IT 


dimension sirj!A(2,2),xn(2),\oc(2),ii?(n) t con(v).v(2,2), 

1 V2(2,2) ,V3(2,2) ,V4(2,2),U1 (2),'J2(2),X(2),VV24(2,2),VVV(2 

2 RH9(4,4),U(4),S(4,4) 

READ (5,1 ) (U(I),I=1 ,4) 

1 FOPriAT (4F1 0.4) 

DO 2 1=1,4 
L=1 


,2),5D(4), 


IF (I.LE.3) L=I+1 

2 READ (5, 3) r.D(I),(RI!0(I,.l),,l=L,4) 

3 F0RMAT(4F10.4} 

DO 4 1=1,4 
L=I+1 


S(I,I)=SD(I)**2 
IF (L.EH.5) GO TO 4 
DO 4 J=L,4 

S(I,J)=Rlin(I,J)*SD(I)*FD(I) 

5(J,I)=S(I,J) 

4 continue 

WRITE(6,5) (U (I) , 1=1 ,4) 

5 FORMAT('T///' MEAT !S VECTOR = ' ,4(F10.4,2X)) 
!!RITE(C,6) 

C FORMAT (' VARIANCE - COVARIANCE MATRIX'/) 

DO 7 1=1,4 

7 WRITE(6,8) (S(I,J),J=1,4) 

8 F0RMAT(5X,4(F10.4,4X)) 

DO 9 1=1,2 

DO 9 J=1 ,2 
VI (I,0)=5(I,J) 

V2( I ,J)=S(I ,J+2) 

V3(I,J)=S(I+2,J) 

9 V4(I,J)=S(I+2,J+2) 

D=V4(1 ,1 )*V4(2,2)-V4(1 ,2)*V4(2,1 ) 

C=V4(1,1) 

V4(l,2)=-V4(l,2)/D 
V4(2,l)=-V4(2,l)/D 
V4(l ,1 )=C/D 
DO 10 1=1.2 
U1(I)=U(I) 

10 U2(I)=U(I) 


L 
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APPENDIX II ( cnruril.TI) ) 


uo n 1=1,2 
do n J=1 ,2 

VV24(I ,J)=0 
DO 11 K=1 ,2 

1 1 VV24 ( I , J )=VV24 ( I , J )+V2( I , K)*V4 ( K,J ) 
DO 12 1=1,2 

DO 12 1=1,2 
VVV(I,J)=0 
DO 12 K=1 ,2 

12 VVV(I,J)=VVV(I,J)+VV24(I,K)*V3(K,J) 
DO 13 1=1,2 

DO 13 J=1 ,2 


13 sinriAd , j)=vi (i , j)-vvv(i ,j) 

COR(l ,1 )=SORT(SinriA(l ,1 )) 
COR(2,2)=SORT(SIfiMA(2,2)) 
C0R(2,l)=Sir;!A(1.2)/(COn(l,l)*COp(2,2)) 
COR(l,2)=COR(2,l) 

WRITE (6 ,1 4 ) SIO HA 
URITE(G,15) COR 

14 FOR’ 3AT ( ' OCOND. VAR. COV .MATRIX • //2 (HX ,F1 0.5)/) 

1 5 FORf 1AT ( ' OSDflCORR . MATR I X ' / / 2 ( 2X , FI 0 . 5 ) / ) 
READ(5,16) II 

16 FOR! ’AT ( 12 ) 

DO 23 n=l,N 


READ (5,1 7) (X(I), 1=1,2) 

17 F0RMAT(2F10.5) 

DO 13 1=1,2 

18 XU(I)=X(I)-U2(I) 

DO 19 1=1,2 
VX(I)=0 


DO 19 J=1 ,2 

19 VX ( I )=VX ( I )+VV24 ( I , J )*XU ( J ) 
DO 20 1=1,2 

20 US(I)=U1(I)+VX(I) 



WRITE(G,21) (X(I),I=1 ,2) 

WRITE(6,22) (US ( I ) , 1=1 ,2) 

21 FORMATCOVALUES CO!!DITIO!IED ON* ,2 ( 5X , FI 0. 4 )// ) 

22 FORMAT( ' CONDITIONAL MEANS ' ,2 (5X ,F1 0.4 )// ) 

STOP 

END 


TRANSFORMATION OF NON-NORMAL MULTIVARIATE DATA 
TO NEAR-NORMAL 

Summary 

A procedure for transforming non-normal multivariate data 
to near-normal data is presented. The procedure is based upon a 
multivariate generalization of a technique proposed by Box and 
Cox (1964). Several examples of the procedure are included along 
with a documentation of the computor software. 

I. INTRODUCTION 

Investigators are often confronted with the problem of 
analysing multivariate data. Upon investigating the existing 
procedures for analysing this type of data, one soon realizes 
that a majority of the existing techniques are restricted to the 
normal distribution. However, real data often violates this 
normality assumption. Thus the investigator is confronted with 
two possible approaches: 1) determine a non-normal multivariate 
distribution which provides a satisfacory model, 2) determine a 
technique for transforming the non-normal data to near-normal 
data. If the Investigator is mainly interested In modeling the 
multivariate data, then the first approach is probably most ap- 
propriate, however, if the main interests are in making statistical 
inferences or probabilistic forecasts then the second approach 
could prove to be adequate. In this paper, we 


have presented a procedure which addresses this second approach. The 
procedure Is a multivariate generalization of a procedure proposed by Box 
and Cox (1964). They proposed the following univariate transformation 
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i 

f 



for X i 0. 



log(y) for X f 0. 


( 1 ) 


Andrews et. al. (1971) extended this transformation to the bivariate case. 

In their paper, they were able to find approximate maximum likelihood 
estimates' for X, by examining the contures of the likelihood function. In 
this paper, the method of Box and Cox is extended to the multivariate case, 
where the maximum likelihood estimate for X is determined using a numerical 
analysis approach. The procedure is presented in a multivariate analysis of 
variance setting, however, several examples are presented which demonstrate 
the versatility of the technique. 


II . Procedure 


Let y. , 


denote a random sample of n- p - dimensional 
n i 


observations from a population with finite mean and finite covariance 
E|, for 1 = 1,2, ..., m. The problem can be stated as; find 

X = (X.j, Xg* .... X/ such that y|^ Is distributed normally with 


* 
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mean ^ 


for i=l 


where 


Since Y 


where 


, and common covarince £ , where 


vU)« (y <X l ) y ( V,T 


( 2 ) 


UJ 

'ij 


' (A.) 

^ y ijk " 1 ^ /A k 
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for A k 4 0 


for A^ ■ 0 


(3) 


(A) 


J»l,2,...,n i , and k=l,2,...,p. For 0, YJj can be written as 


4*3 - D ' 1(Y iJ - J> 


(U) 


D 5 diag(A^,Ag, . . . ,A^) 

J is a pxl vector of l's 


j A. Ap A rp 

Y iJ " ^ y ijl’ y i32’'**’ y ijp ) ' 


(A) 

ij 


its density function can be written as 
f(z) = exp{-l/2(z-u) T £" 1 (z-u)}(2ii)" p/2 |£r 1/2 (5) 


, = Y^. From this, one can determine the density function for the 
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untransformed data w ■ as g(w) f(z) where. 


az 


IJ ' k ?, “S * Jl (1 V 


V 1 


k*l 


( 6 ) 


Hence the joint likelihood function becomes 

LW-(; n1 n K,j) (2n)' np/Z \l\' n/i 
1*1 j=l 

m n i tt) T , (X) 

* exp { -Js Z Z (Y.. - y) 1 Z (Y>. - u) > 
1=1 j=l iJ 


where n 


m 


Z n, . The likelihood function can be written as 
1=1 1 


(7) 


L(X) = K(2ir)' np/2 |E| _n/2 exp {-np/2> 


? n. 


(8) 


where K = n n K . . and u and Z are replaced by their maximum likelihood 
1=1 j=l 1J 

estimates 

a -j n^ (X) „(X) 

V, = 2 Y.. — 1 Y, 

1 n f j=1 Ij 1 


m n 


Z _ 1 Z- Z (Y.. 

"n1=l J-l ij 


1 (X) .(X) (X) Jx) 

- *! ) (Y^ - Y 1 ) T . 


( 9 ) 
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Equation (8) follows from equation (7) since 


m n. 


1 (X) A T (X) A 

1 i (Y^ - yj) e (Y,| j - Wf) 

(X) (X) 


1-1 j -1 
m n 


1 a , (X) (X) 

» E E (tr E (Y. . - y.) (Y . - U.) T ) 

1-1 j-1 1J 1 1J 1 


. . m (X) 

7 * / ** /u 


(X) 


- tr E ' ( i i'( Y.. - J (Y - y,)' ) 

1=1 j-1 1J 1 U 1 

A * A 

= n tr (E~ E) = np. 

Equation (8) can be further simplified as 

L(X) = C • h(X) (10) 

where C = (2n)” n ^^^ exp (-np/2) 

h(X) = |K~ 2/n E |" n/2 (11) 

Note that maximizing the likelihood function L(X) is equivalent to 
minimizing the function h(X) This function can be further simplified by 
considering 

K 2 '" - S 1 k|< ) 2 ^ 

i=i j-i 

- I ( n X 1 (y ijk ) V1 ) 1/n ) 2 


k=l 1-1 j -1 


P . X.-l 2 

,n (y k ) k ) z 

k=i K 


( 12 ) 
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• m rij i/ii 

where y. ■ ( IT I! y <4t> ) is the geometric mean for the k variate, 
k 1b1 jsl UK 

A 

k * 1, 2, p. From equation (4) E can be written as 


1/n 


(X) (X) (X) (X) T 

2 - 2. £ (Y^ - 7 1 ) (Y<< - 7 ) 


1=1 j=l 


1j 


m n 
= E E 


1 n"l 


X T 


1=1 j=l 


D < Y 1d - Y 1 > < Y ij - V 


■rl 


(13) 


Hence |e|, becomes 


a >, m n. a a a a i 

1 2 1 = | D~ ' | | Z E 1 (Y.. - Y. ) (Y - Y. )| |D"'| 

1=1 j=l 1J 1 1J 1 


* |0' 2 | |s| 


(14) 


m n, X \ \ T 

where G = E E (Y^ - Y^) (Y^ - Y^). Thus the minimization of 


1=1 j=l 


-1 


h(x) Is equivalent to minimizing 


0(X) = 


M 


j k 2/N d 2| 


M 


p V 1 2 
( n x.y. k ') • 

k=l k K 


(15) 


n X X 


Note that equation (15) reduces to E (y. - y Y 

1=1 1 


X-l 2 

(tf ) 


( 16 ) 
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which was proposed by Box and Cox (1964) for the univariate case. 

The function 0(X) In equation (15) can now be minimized using a 
standard numerical technique. In this paper the Flecher-Pcwell algorithm 
of deflected steepest descent Is used (see Appendix A). 

III. Application 

The first example Illustrates a violation of the equality of covariance 
matrix assumption in a multivariate analysis of variance problem. The data set 
is R.A. Fisher's classical iris data (Fisher, 1936) where the response 
measurements are sepal length, width and petal length, width for three iris 
species: virginica, versicolor, and setosa. Although this data was originally 
presented as an application of linear discriminate analysis, Morrison (1967) 
uses this as an example in multivariate analysis of variance, for which he 
states, "we shall of course assume... a common covariance matrix". However, 
in applying Bartlett's likelihood ratio test for equality of covariance, we 
obtain a test statistic of 141 for 20 degrees of freedom. Hence the hypothesis 
of equality of covariance can easily be rejected with a high level of 
significance. In figure 1, the confidence ellipse for the two untransformed 
variables: sepal length and sepal width, clearly illustrate the difference 
in covariance matrices. The data is then transformed, and the corresponding 
confidence ellipses are presented in figure 2. Although the confidence 
ellipses for the transformed data are more nearly identical, Bartlett's test 
statistic has been reduced to 63, however, this value is still significant 
at the .01 level. 


Untransformed 95% Confidence Ellipses 
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Sepal 

Width 


Figure 2. Transformed 95% Confidence Ellipses 
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In the second example, we are Interested In obtaining probabilistic 
forecasts. The data was originally preserved In a paper by Haggard et. al. 
(1973), In which the author was able to model the maximum rainfall from 
tropical cyclone systems across the Appalachians using the Gamma distribution. 
Since or. a of their primary objectes was to obtain estimates for the 
probability of rainfall exceedence In the Appalachian regions, I felt that 
comparative results could be obtained by transforming the data then using 
the well tabulated normal distribution. The results are given In Table 1. 

IV. Conclusions 

A method transforming non-normal multivariate data to nearly-normal 
data Is presented. The method extends the univariate transformation of 
Box and Cox (1964) to the multivariate case. A numerical method for 
approximating the optimal transformation Is also Included (see Appendix A). 

The procedure was then applied in two applications. The first was In the 
area of multivariate analysis of variance where the primary objective was to 
achieve equality of covariance matrices. It was shown that the transformed 
data was less heterogeneous than the un transformed data. However, the 
population covariances were still unequal. The second application Illustrated 
that this type of procedure can be used when the primary objective Is the 
estimation of tall probabilities. This method allows the use of the normal 
distribution on the transformed data, rather than determining the appropriate 
non-normal distribution for the untransformed data. 
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TABLE 1 

Expected Probabilities of Exceeding Arbitrary 
Precipitation Amounts Over the Appalachian Region 


Precipitation Data Set* 

In Inches 

ABC D 

Hi 



I 

II 

I 

II 

I 

II 

I 

II 

1 

.978 

.966 

.993 

.999 

.981 

.971 

.995 

.997 

2 

.913 

.903 

.959 

.999 

.932 

.924 

.971 

.976 

3 

.821 

.819 

.893 

.962 

.865 

.864 

.926 

.931 

4 

.717 

.723 

.806 

.809 

.789 

.794 

.866 

.866 

5 

.613 

.624 

.706 

.625 

.710 

.719 

.794 

.788 

6 

.515 

.528 

.605 

.472 

.631 

.644 

.717 

.706 

7 

.427 

.439 

.507 

.361 

.556 

.571 

.639 

.623 

8 

.349 

.359 

.418 

.283 

.486 

.500 

.562 

.544 

9 

.283 

.291 

.340 

.227 

.423 

.436 

.489 

.471 

10 

.227 

.232 

.273 

.186 

,365 

.376 

.422 

.405 

15 

.070 

.066 

.079 

.090 

.165 

.166 

.182 

.174 

20 

.019 

.016 

.020 

.057 

.070 

.066 

.070 

.076 

25 

.005 

.003 

.002 

.042 

.028 

.025 

.025 

.032 

30 

.001 

.001 

.001 

.033 

.011 

.009 

.008 

.023 


* A - maximum 24-hour precipitation all storms. B - maximum 24-hour 
precipitation from no more than one storm per year. C - maximum 
precipitation totals from all storms. D - maximum precipitation totals 
from no more then one storm per year 

«« 

I- gamma parameters from Haggard et.al. (1973); II transformed 
normal probabilities. 
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Appendix A 

Flecher-Powell method of deflected steepest descent, requires the 
gradient vector 


0(x) = 


where 


_30 

3X 

30 

3X 


30 

3X_ 


(A.l) 


0(X) = 


M 


P « X^-l 
( n X.y. ,2 

k=l k k ; 


(A. 2) 


3 0(X ) 
3X, 


h 8A h k=l 


P X. -1 0 

D X -1 8(11 Vk k } " 2 l 6 l 

< " Vk > + - JL± TT (A. 3) 


/ p • V \-2 

3( n x k y k ) 


k=l 


3X. 


-2( n X y. k ) 2 
k=l K K 


fX t, + ln V X n 


-1 


(A. 4) 


Since |G| = l g. s a.- -• where 
j=l 1J 


G = (fij) 


(A. 5) 


a^j is the cofactor of g^j 
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Also, since g^j only depends upon X^, Xj using the chain rule we have 


liilL 

3*u 


. I J |1£L 

1=1 j-1 3g ij 3A h 


(A.6) 


where 


ii£L 

39 1j 


Hj 


(A. 7) 


and 


99 


uv „ 


3X. 


< b 2 


and 


If u, v t h 
If u or v = h 
If u = v = h 


a m n, X n X n X. X. 

b 2 = 3X^~ ^ (y iju ' y 1u } (y ijh ‘ y 1h* ^ 


(A. 8) 


m n^ X X X. 

= i 3^(y. u ■■ u ' h 


l l'(y. u . - y. u ) (y.'l.ln y. .. - y b lr, y. .. ) 
,j = 1 j-i iju iu ’ VJ ijh ■ 7 ijh •'ih ■ 7 ijh' 


id n « Xl 


b, = 2 1 E 1 (y. 


— A h 


1=1 j=l 


1 jh 


y 1h > (y i jh 1 " y ijh- y i!> y 1Jh> 


(A.9) 


From this, equation (A. 3) becomes 


90(X) _ 


P * A k ~^ \2 
( n x.y. k )* 

k=l K K 


P 

k=l 

m 


a, 


kh 


9g kh + fchh !5]h 


3X, 
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Test of Fit for the Extreme Value Distribution 
Based Upon the Generalized Minimum Chi-Square 

Summary 

A goodness of fit test for the extreme value distribution is developed. 
The procedure is based upon the generalized minimum chi-square distribution 
[Gurland and Dahiya (1579)] . Application of the test is given for some 

extreme value data [Gumbel (1964)]. 

I. Introduction 

There are several difficulties with using the Pearson chi-square test 
of fit for continuous distributions [c. f. Dahiya and Gurland (1970 )]. 

These difficulties are primarily concerned with the choice of cell width and 
the number of cells. However, to the applied statistician or non-statistician 
who must use test of fit procedures on a frequent basis, the primary difficulty 
of the procedures is in the users set up. That is, the user must have knowl- 
edge of the tabular values for the null hypothesis. Dahiya and Gurland 
(1970 , 1972) presented a goodness of fit test for several continuous dis- 
tributions which eliminate most of the user’s set up. Their procedure was 
based upon the generalized minimum chi-square statistic. In this paper, I have 
developed a test of fit for the extreme value distribution based upon this 
generalized minimum chi-square technique. 
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II. Procedure 

Suppose that one vould like to test the null hypothesis given by 

H q : X lt X 2 X n * F x (x;8) (1) 

where X^,Xg,...,X n denotes a random sample of n observations from a distribution 
function F x (x;0). F^ is an asymptotic Fisher-Tippett type 1 distribution, that is. 


F x (x;6) = exp{ -exp( -( x-a)/0 ) } (2) 

_«o < a < “ 

6 > 0 . 


Let T denote a transformation from the population raw moments to 5, which 
can be written as a linear function of the parameters 6 where 


T 

n' = (nj^rig, . . . ,Hg) 

€ = U ± ,i 2 5 s ) T 


(3) 


+■ Vi 

and nl is the j raw population moment for F Y and £ = W0 , V is a known sx2 

J fp ^ 

matrix, and 0 = (a, 6) . That is, 

T: n £ = W8. (M 

Let m = (ra 1 ,m 0 ,...,m ) denote the sample raw moments corresponding to n 

Jl tz S 

T 

and let h = (h 1 ,h_,...,h ) denote the sample values corresponding co £, 

JL c S 

f 

T: m -*■ h. 


that is, 


(5) 


By the central limit theorem, we know 
n(m' - n * ) ^ n(0, G) 

where the 1j^ element of the matrix G is 


g.-. = n' - n’ n* 
1J i + j i j 


for 1 , j = 1 , 2,. . s. It also follows that 
n(h -s) n, n(0 ,E) 


where E = TGT . Now using the distributional properties for the quadratic 
forms, we know that 

Q* = n(h -S) T E" 1 (h -S) (9) 

has an asymptotic chi-square distribution with s degrees of freedom where 

/\ 

E is a consistent estimator for E. Since S = W0, an estimate for 6 can be 
found by minimizing Q*. In which case, the estimate becomes 

6= (W T E -1 W) -1 W T r 1 h. (10) 

/s /\ 

By letting S = W9, Q* becomes 


where 


Q = nh Ah 


A * E~' (I - R) 


^ T A -1 _ 1 T A - 1 

R = w(w' l 'w) V E 


Again by the distributional properties of the quadratic forms, 0 has a non- 

aa 

central chi-square distribution with degrees of freedom = tr EA and non 

y A A 

centrality parameter X = S AS if and only if EA is an idempotent matrix. 

AA A AA ^ 

It is easy to verify that (EA)^ = EA, and X = 0, so Q has a chi-squr.re 


. 4 * 
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distribution with s-q degrees of freedom. Using this distribution, one 

can reject the null hypothesis (1) with type I error if Q > *( s-q) , where 

x o 


Pr( X >_ x^(s-q)) ■ a. 


(13) 


Dahiya and Gurland (1970) developed the non-null distribution for Q, using this 
distribution one can compute the power of the test for a specified non-null 
distribution. In order to test (l), the transformation T and the matrix W need to 
be specified. Since we know that the populations cumulants for the extreme value 
distribution are 

•c, ■ (-B)V J ~ 1) for J =2,3,... (lU) 

J (1) 


where 


= (-l) n+1 n!6(n+l) 

( 1 ) 

6 ( n ) = “ i -n . 

i=l 


(15) 


-1 -1 -1,T . .. / (2) . (l) ,(s+l), ,(s)vT 

By letting £= (k,k 0 »•••»< ... 9**1 ) and W= (^ /<J> /* ) , 

3 2 4 3 s+2 s+1 (1) (D (l ) (1 ) 

I 

and 6 = B it is possible to map n •* £ where s = 4 and q = 1. By letting 

T i 

h = (h^jhgjh^jh^) , where hj = k j + 2/ k j+i» for J=l,2,.., k , and is the J 

sample cumulant. We are now able to compute Q, where 


.th 


z = JGJ T | 6= ~ (16) 

H 

J = (j ); J = T~° for m,n=l,2, . . . ,s 
mn mn die 


and B is the maximum likelihood estimate for 6 . 


The values In equation (15) can be found In Abrahomovich, hence J becomes 


1/6 


10 0 0 
-.885 .6079 0 0 

0 -1.131 .4174 0 

0 0 -.5901 .154 


(17) 


From these values, we are able to compute Q in (11) for the sample values 

* ? 

X.j, X 2 , .... X n . Hypothesis (1) can be rejected if Q > x (3) since 
s * 4, q « 1. 


Applicati on 

In this section, an extreme value data set given in Gumbel and 
Goldstein (1964) is analysed using this test of fit procedure. The data set 
consists of the oldest ages at death for men and women in Sweden from the period 
1905-1958. The data for male and female are fitted separately. Gurrbel and 
Goldstein (1964) estimated the extreme value distribution parameters using 
a modified method of moments. Tables 1 & 2 contain a comparison of the two 
different procedures in term of estimated parameters and cumulative tail 
probabilities. It must be noted, that the null hypothesis of the extreme 
value distribution being the null distribution could not be rejected at a 
significance level of greater than 70%. 

In th* second example, extreme monthly temperatures and winds for 
three United States locations were analysed. The data set taken from th« 
dally meteorological records, 1970-1971, for New Orleans, LA., Orlando, 

FL., and Daytona Beach, FL. The results are summarized in Tables 3 and 4. 
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Table 1: Comparison of Procedures using Swedish Men 


Method of 2 

Moments Generalized minimum x 


A 

a 

A 

e 

X* 

F x (x) 

a 

a 


to > 

X 

* 

G x (x) 

102.49 

1.39 

100.90 

.0433 

102.53 


1.25 100.90 

.0251 



101. 6C 

.1625 



101.66 

.1346 



102.61 

.3994 



102.61 

.3914 



103.24 

.5582 



103.24 

.5674 



104.22 

.749 7 



104.22 

.7720 



105.72 

.9067 



105.72 

.9250 



106.50 

.9457 



106.50 

.9591 


* the values X 

represent the 5, 

10, 20, 

30, 

40, 50, 54th 



smallest sample 

value. F x and 

G x are 

the 

corresponding c.d.f. 




Table 2: 

Comparison 

of Procedures using Swedish Women 



Method of 



2 



Moments 



Generalized minimum x 


A 

a 

A 

6 

X* 

F x (x) 

A A 

a 8 X' 

G x (x) 

103.83 

1.25 

102.54 

.0604 

103.33 1.57 102.54 

.2118 



103.31 

.2196 

103.31 

.3866 



103.94 

.4002 

103.94 

.5293 



104.52 

.5623 

104.52 

.6442 



106.15 

.8553 

106.15 

.8558 



106.50 

.8889 

106.50 

.8829 


* same as in Table 1 
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TABLE 3 

Extreme Monthly Temperatures 


Site Extreme Value Distribution 



a 

B 

Q* 

New Orleans 

83.8 

.98 

.001 

Orlando 

84.8 

.88 

.003 

Daytona Beach 

81.7 

.67 

.002 


* null distribution of extreme valued distribution can not be rejected. 


TABLE 4 

Extreme Monthly Winds 


Site Extreme Value Distribution 



A 

a 

A 

6 

A. 

Q 

New Orleans 

15.4 

2.9 

.9 

Orlando 

13.6 

?.’S 

.6 

Daytona Beach 

13.0 

2.2 

.,4 


* same as In Table 3 


IV. Conclusions 


A procedure for testing the goodness of fit for the extreme value 
distribution, based upon a generalized minimum chi-square is presented. The 
procedure is applied to several data sets where the extreme value distribution 
is a potential fit, although it must be mentioned that the meteorological data 
set was included in a manner which lends itself to program utility rather than 
for meteorogical interpretation. 
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Test of Fit for Continuoas Distributions Based Upon 
the Generalized Minimum Chi-Square 


Summary 

A procedure for test of fit for several continuous probability 
distributions based upon the generalized minimum chi-sqare method is 
presented. The procedure was first presented in a series of papers by 
Dahiya and Gurland ( (1970a) , (1970b) , (1972) ). Examples of the procedure 
are included, along with the corresponding computer listing. 

I Introduction 

Dahiya and Gurland (1970a) discuss the difficulties with using the 
Pearson chi-square test of fit for continuous distributions. These dif- 
ficulties are primarily concerned with the choice of cell numbers and widths. 
However, to the applied statistician who must use test of fit procedures on 
a frequent basis the main disadvantage is in the users setup. That is, the 
user must have knowledge of the parameters and the tabular values for the 
specified null distribution. These demands severally hamper the investigator 
who must determine an appropriate distribution from potentially many distribution 
functions. The purpose of this paper is to present a test of fit for continuous 
distributions wi.ich minimizes the users interface in the estimation of 
parameters for the specified null distribution or in specifying the tabular 
values of the null distribution. In fact, several different families of 
distributions can be tested for fit using a single 
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setup. The procedure Is based upon the generalized minimum chi-square 
(GMCS) statistical method. Section 3 contains the GMCS procedure for the 
univariate normal and gamma distributions. 


Procedure 

Suppose that we want to test the null hypothesis 

Ho: x 1 , x 2 , . . .x n -V F x (x;e) « J(x;o) (1) 

where x^, Xg,..^ Is a random sample of n-observatlons from an unknown 

distribution function F x (x;e); a Is a q x 1 vector of parameters and 

l?(x;0) Is a specified family of distributions with admissable parameters 0. 

The (GMCS) procedure can be used for testing any family of dis- 
tribution ;f(x; 0 ), provided there exists a transformation T, where 

T: y - C (2) 

where u' = (y * ,y' ,. . .u ' ) T , y' Is the raw population moment 
12 s j 

and £ = (^l * * * * ) T can * 5e ex P ressed as £ « We (3) 

fora known s x q matrix w and s > q. Let m' * (m' , m',...m') T 

1 2 s 

denote a s x 1 vector of raw sample moments and define 
h B (h^, hgi.-.h^ to be the image of the transformation T, that Is 

T: m* •+ h. Using the central limit theorem, we have 

n(m' - y') - n(0, G) (4) 

where G * lg<J, g, < * y‘ - y' y', 1, j ■ 1, 2, ...,s. 

1J 1 + j i j 
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From this. It can be shown that 

n(h - £) -N(P, E) (5) 

where E ■ ii(5j\ J the Jacobian matrix for the transformation T. Now using 
the properties of quadratic forms, we know that 

Q - n(h - 0 T £ -1 (h - 0 (6) 

has a chi-square asymptotic null distribution with s degrees of freedom. 
Furthermore, t.hls distribution does not change when we estimate E In (6) 

A A 

by E, where E Is a consistent estimator for E. Since £ « We, we can 

A 

estlma ;e 6, by finding fl which minimizes Q. This estimate Is given by 

0 - (W T E^W)* 1 W T E _1 h. (7) 

By letting £ * Me, the minimal Q Is 

Q ■ n(h» 0^ "\h - £) s nh^Ah (8) 

where 

A A 1 A 

A « E“'(I - R) (9) 

R « W(W' E 1 W) 1 W 1 . 

A 

Again, using the properties of the quadratic forms, we know that Q has a 

A A 

non-central chi-square distribution with degrees of freedom ■ tr(E A) and 

TA A A 

asymptotic non-centrality parameter X « £ A£, If and only If EA Is Idem- 

A A 

potent. Under the null hypothesis, tr(EA) * s - q and X ■ 0. Hence the 

A. A 

asymptotic distribution of Q Is x (s - q). Using this distribution, we 

A 

can reject the null hypothesis with a type I error If Q > x a (s - q). 

A 

Gurland and Pahlya (1970) developed the non-null distribution for Q. 
Using this result, they were able to compute the power of the test for 
selective alternative distributions. 
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1 0 
C 1 
0 0 
0 2 

The transforamtion T from n to ( can be achieved in two steps ; : p p 

Tg: p Hence, E in equation (5) becomes 

E « J 2 J 1 CJ 1 J 2 (13) 




. 1 - 


o » ® 
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where 

V 

Ji * (j m ); j * --j— — m, n « 1, 2 s 

M 

e 5 ■ % u, v = 1, 2, .... s. 

* uv 

By assuming that u' * 0, J, and J, become 

1 1 * 

» 

1 0 0 

0 1 0 

°1 = -30 2 0 1 

0 0 0 

1 0 0 

0 l/0 2 0 

0 0 1 

0 0 0 

and equation (14), becomes 

e 2 0 0 0 

0 2 0 4 

0 0 66 * 0 

040 32/3 






05) 


( 17 ) 
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Since 6g Is unknown, let 0 2 denote the usual maximum likelihood estimate. 


Then l s E 


Now by computing 0 and Q In equation (7) and (8), 


e 2 ■ * 2 . 


one can test the hypothesis (10). 


Gamma Distribution 


Test the hypothesis 

Ho: X 1# Xg,..., X n *F x (x;0) * r( 0 ] , 0 2 ) 
where the density function for the gamma distribution r(e.j, 0 2 ) is 

0,-1 


f y(x » 0^ ,©2 ) - 


e " y y 1 
0 2 r(e 1 ) 


; y = x/© 2 


(18) 


(19) 


©■j ,0 2 > 0. 

Since ^ = (j - 1)! 0^0^, the cumulant, we can express £ = W0*, 


where 


£ = ( -j »? 2 - ^ 2 ~ < 3<2 • £4 - k 4 k 3 ) 


(20) 


W = 


1 

0 

0 

0 


0 

1 

2 

3 


0 * - (o* = 0 ^ 2 .o 2 = e 2 ) 


( 21 ) 


The transformation T from n' to 5 can be obtained in two steps 


T ] : n' + K 


T 2 : K 


(22) 
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where < 

“ (< 

1 » k 2 ’*3 

.< 4 ) T * 

In which case Z 

becomes 

Z 

* J 


V 






where 





* 






1 

0 

0 

0 






, % 

1 

0 

0 





s 

% 

j 23 

1 

0 






j 14 

J 2 4 

j 34 

1 





j 12 

= -2 

"i 




J 13* 

-3n ' 
1 


j 23 

= -3n 

2 + 6n l 




in’ 

-4n 


J 2 4 

» -6 

' z + 12(np 2 



d 3 4 * 

S 


n i 

= r( 0 l 

+ j) 

j 

0 2 * 

3 s 

1, 2, 

3, 4, 

• • • 


(23) 


(24) 


. n 3 


*4^3 + ^^3 “ 24(Tvj ) 


(25) 


r(0j 


-Ko< 


-1 -1 


2*1 


0 

0 


-1 

‘*3*2 


~V3 


0 
0 
0 

■1 .-1 


(26) 
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Since e i»02 are unknown, they can be estimated by e-| ^ where 

e 2 - x /© 1 

e 1 - y" 1 / 4 (1 + (1 + 4y/3) % ) 
y s log (I/GM) 

7_1 ? X i (27) 

x " n 1=1 1 

n 1/n 
GM = (n X.) 

1=1 1 

A A A 

By replacing In E, we can test the hypothesis (18) using Q. 

Results 

In order to demonstrate the GMCS procedure, the procedure was used 
in three different experiments. The first was to simulate data from several 
different distributions and determine the test of fit. In the second example 
the procedure was analysed using meteorological data consisting of several different 
atmospheric variables. The third experiment consisted of analyzing a meteorological 
data set from a specified distribution function. 
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E xperiment 1 

In this experiment, random observations were simulated from many 
different distribution functions In order to demonstrate how robust the 
procedure Is to varylnq sample sizes, shape parameters, etc. This part of 
the experiment was not meant to provide conclusive evidence that the 
(6MCS) procedure Is better or worse than any other procedure, but was 
Intended to point out any apparent deficiencies. The results have been 
summarized In Table 1. In this table, I have only Included the results 
for fitting the true distribution, however, the procedure may have Indi- 
cated that another distribution could have provided satisfactory fit. 
However, this Is explainable since the Gamma and Extreme Value distribution 
can resemble many other distributions depending upon their shape parameters. 
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TABLE 1 

Evaluation GMCS procedures using Simulated Data 


True 



Sample 

Estimated 

A 

Distribution 

Parameters 

Size 

Parameters 

Q 

T(Y t P) 

y 

$ 


A 

y 

A 

6 



3 

1 

10 

1.2 

1.06 

6.600* 


II 

II 

25 

1.0 

.91 

.001 


It 

II 

50 

1.1 

.86 

1.200 


II 

II 

100 

1.1 

.90 

4.900* 


2 

1 

10 

.97 

.98 

5.600* 


II 

II 

85 

.88 

.88 

14.900* 


II 

II 

50 

1.18 

.96 

12.700 


II 

II 

100 

.83 

.72 

3.300 


.5 

1 

10 

1.99 

1.6 

.420 


<1 

II 

25 

.80 

.63 

1.000 


11 

II 

50 

1.02 

.77 

1.200 


II 

it 

100 

1.17 

• 91 

15.100* 

N(y, a 2 ) 

V 

2 

a 

N0B 

A 

y 

"2 

a 

A 

Q 


10 

25 

10 

12.1 

11.8 

.008 




25 

9.5 

31.5 

.091 




50 

8.9 

20.0 

.041 




100 

10.2 

23.9 

.001 

Extreme 

a 

6 

NOB 

A 

a 

A 

B 

A 

Q 

value, a, 3 

5. 

1. 

10 

5.01 

1.68 

.001 




25 

5.04 

1.15 

.008 




50 

5.04 

.85 

.003 




100 

4.82 

.85 

.006 


2. 

2. 

10 

2.90 

.98 

.002 




25 

2.69 

1.50 

.004 




50 

1.74 

2.08 

.017 




100 

2.09 

1.95 

.033 

Exponential 


X 

NOB 

X 


A 

Q 

X 


.5 

10 

.69 


9.04 * 




25 

.56 


3.2 




50 

.53 


1.3 




100 

.42 


4.15 * 



1.0 

10 

1.04 


.33 




25 

.83 


.75 




50 

1.28 


.29 




100 

1.1 


2.90 



2.0 

10 

2.60 


4.9 * 




25 

1.89 


1.54 




50 

1.97 


1.04 




100 

1.92 


.35 


* null hypothesis can be rejected at a e 0.5 level 
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M 




In this experiment meteorological data sets from three southern 
United States locations were analysed. The first set consisted of monthly 
percepltatlon totals and monthly mean tenperature for the years 1936-1975 
for sites New Orleans, LA, Orlando, FL, and Daytona Beach, FL. The results 
for these data sets have been summarized In Tables 2 & 3, where the data 
sets are partitioned Into five year Intervals, each containing 60 
observations. The second data set consists of dally (high temperature, 
maximum wind speed) for the three U.S. sites. The observations are 
partitioned Into monthly intervals for the 1970-1971 data. The results 
are summarized In Tables 4 A 5. Tables 6 & 7 contain the results for test 
of fit for extreme monthly temperature and wind for the three U.S. locations. 

It should be mentioned that the above data set was partitioned for the 
author's convenience rather than for meteorological interpretation. 
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TABLE 2 

Monthly Total Precipitation 


Site** 

Year 

A 

V 

Normal 

7 

A 

Q 

A 

X 

Exp 

A 

Q 

A 

Y 

Gamma 

A 

A 

Q 

Extreme 

A A 

a S 

I 

1936-40 

4.8 

24 

3.3 

.02 

7.6* 

1.9 

.04 

.0 

12.9 

3.0 


41-45 

4.5 

12 

1.3 

U 

8.2* 

2.2 

.05 

II 

8.3 

2.4 


46-50 

5.4 

20 

.8 

II 

5.3* 

1.8 

.03 

II 

9.4 

3.1 


51-55 

4.6 

9 

1.4 

II 

7.8* 

.... 



6.7 

2.2 


56-60 

4.7 

8 

.3 

II 

10.6* 

2.7 

.06 

II 

7.3 

2.1 


61-65 

4.6 

7 

.0 

II 

9.3* 

2.3 

.05 

II 

6.3 

2.1 


66-70 

4.4 

8 

.3 

II 

8.6* 

2.3 

.05 

II 

6.7 

2.2 


71-75 

5.8 

10 

.3 

II 

13.1* 

3.5 

.06 

II 

8.8 

2.4 

II 

1936-40 

4.2 

13 

.7 

.02 

3.7* 

1.6 

.04 

.0 

7.5 

2.5 


41-45 

4.0 

14 

1.4 

II 

3.6* 

1.6 

.04 

1) 

8.8 

2.4 


46-50 

4.5 

19 

.5 

ii 

.9 

1.1 

.02 

II 

7.7 

3.0 


51-55 

4.4 

16 

.9 

it 

.9 

1.5 

.04 

II 

8.0 

2.7 


56-60 

3.4 

9 

.1 

ii 

2.0 

.... 



5.3 

2.2 


61-65 

4.3 

19 

1.2 

n 

1.4 

1.3 

.03 

It 

9.0 

2.9 


66-70 

4.0 

9 

1.4 

ii 

4.6* 

1.7 

.04 

II 

6.1 

2.2 


71-75 

3.9 

15 

1.2 

M 

1.3 

1.2 

.03 

II 

7.8 

2.6 

III 

1936-40 

3.8 

7 

1.5 

.02 

6.1* 

1.8 

.05 

.0 

5.6 

2.0 


41-45 

4.5 

16 

.3 

11 

2.0 

1.3 

.03 

II 

7.4 

2.8 


46-50 

4.4 

13 

.3 

II 

3.5* 

1.5 

.03 

II 

7.1 

2.6 


51-55 

4.1 

20 

1.9 

It 

2.2 

1.3 

.03 

II 

9.3 

2.8 


56-60 

3.9 

10 

.3 

II 

3.3* 

1.5 

.03 

II 

6.2 

2.3 


61-65 

3.9 

9 

.3 

II 

9.9* 

1.7 

.04 

II 

6.1 

2.2 


66-70 

3.9 

14 

.8 

II 

1.3 

1.2 

.03 

II 

7.3 

2.6 


71-75 

3.9 

9 

.3 

II 

5.3* 

1.7 

.05 

II 

6.0 

2.2 


* null hypothesis can be rejected at oc = .05 level 

** I - New Orleans; II - Orlando; III - Daytona Beach 


Q 


3.5* 

1.8 

3.3* 

1.4 

1.1 

1.1 

1.2 

1.3 


2.2 

2.3 

3.9* 

2.7 
1.9 
3.5* 
1.5 

2.8 


1.2 

3.0 

2.3 

3.5* 

1.9 

1.5 

2.7 

1.5 
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TABLE 3 

Monthly Mean Temperature 


Site** 

Year 

A 

V 

Normal 

; 2 

A 

Q 

A 

X 

Exp 

A 

Q 

A 

Y 

Gamma 

A 

6 

A 

Q 

A 

a 

Extreme 

A 

6 

A 

Q 

I 

1936-40 

69.7 

116 

.0 

.01 

24.1* 

26 

.36 

.0 

75.2 

8 

.0 


41-45 

69.4 

120 

it 

II 

II 

25 

.39 

H 

75.1 

9 

ii 


46-50 

69.4 

106 

it 

II 

II 

29 

.39 

II 

74.6 

8 

ii 


51-55 

69.2 

111 

it 

It 

II 

25 

.42 

II 

74.9 

II 

it 


56-60 

68.5 

121 

II 

II 

II 

25 

.36 

It 

74.2 

II 

ii 


61-65 

67.5 

121 

H 

II 

II 

22 

.32 

II 

73.0 

It 

ii 


66-70 

67.0 

130 

H 

II 

II 

23 

.34 

II 

72.9 

9 

n 


71-75 

68.6 

99 

it 

II 

II 

23 

.49 

II 

73.8 

8 

n 

II 

1936-40 

71.0 

69 

.0 

.01 

24.6* 

40 

.57 

.0 

75.2 

6 

.0 


41-45 

72.0 

80 

It 

II 

II 

41 

It 

II 

76.7 

7 

It 


46-50 

73.4 

58 

II 

1* 

II 

43 

II 

It 

77.0 

6 

II 


51-55 

71.8 

72 

II 

II 

II 

57 

.80 

II 

76.0 

7 

II 


56-60 

71.8 

78 

It 

II 

II 

34 

.48 

II 

76.1 

7 

It 


61-65 

72.4 

73 

II 

11 

It 

40 

.55 

II 

76.6 

11 

II 


66-70 

71.8 

83 

II 

II 

II 

36 

.51 

II 

76.3 

II 

II 


71-75 

73.6 


II 

II 

II 

53 

.72 

It 

77.4 

6 

II 


III 

1936-40 

69.7 

63 

.0 

.01 

24.7* 

41 

.54 

.0 

73.7 

6 

.0 


41-45 

70.1 

86 

Ii 

II 

II 

33 

.47 

II 

74.8 

7 

II 


46-50 

71.5 

61 

II 

II 

II 

40 

.56 

II 

75.3 

6 

II 


51-55 

70.4 

75 

11 

II 

II 

55 

.78 

II 

74.9 

7 

II 


56-60 

70.0 

82 

II 

II 

II 

32 

.46 

II 

74.5 

7 

II 


61-65 

69.8 

76 

II 

II 

II 

39 

.56 

II 

74.3 

7 

It 


66-70 

70.0 

89 

11 

II 

li 

34 

.49 

II 

74.7 

7 

II 


71-75 

71.3 

60 

II 

II 

it 

34 

.50 

II 

75.2 

7 

It 


* null hypothesis can be rejected at a » .05 level 

** I - New Orleans; II - Orlando; III - Daytona Beach 


100 


TABLE 4 

Dally Maximum Temperature 


. SI te** Date*** 

Normal 


Exp 



Gamma 



Extreme 



.a 

y 

A 2 

a 

A 

Q 

A 

X 

A 

Q 

a 

y 

a 

6 

A 

Q 

A 

a 

A 

6 

A 

Q 

I 1/70 

47.5 

141 

.0 

.02 

11.6* 

18.9 

.39 

.0 

55.8 

9. 

.0 

3/70 

60.0 

44 

m 

.01 

12.8* 

49.2 

.73 

ii 

63.8 

5. 

.0 

6/70 

79.7 

20 

it 

II 

II 

116.0 

1.4 

ii 

81.8 

3. 

.0 

10/70 

69.0 

26 

it 

II 

II 

99.8 

1.4 

it 

71.7 

4. 

.1 

1/71 

55.3 

112 

ii 

It 

II 

23.5 

.42 

it 

61.0 

8. 

.0 

3/71 

59.4 

75 

ii 

II 

II 

23.3 

.38 

ii 

64.6 

7. 

.0 

6/71 

80.2 

5 

ii 

II 

II 



.... 

81.9 

2. 

.0 

10/71 

71.8 

23 

ii 

II 

II 

70.5 

oc 

O'. 

• 

H 

74.1 

4. 

.0 

II 1/70 

55. 

94 

.0 

.01 

12.2* 

18 

.32 

.0 

60. 

7. 

.0 

3/70 

76. 

22 

11 

II 

II 

11 

1.5 

It 

78.4 

4. 

11 

6/70 

83.9 

1 

II 

II 

II 

67.8 

8. 

II 

84.4 

1. 

II 

10/70 

63.5 

66 

II 

II 

II 

22.3 

.35 

II 

67.1 

6. 

II 

1/71 

64.3 

87 

II 

II 

II 

22.2 

.34 

II 

68.8 

7. 

II 

3/71 

72.0 

53 

II 

II 

II 

54.1 

.74 

II 

76.0 

6. 

II 

6/71 

83.4 

3 

II 

II 

II 

45.8 

5.5 

It 

84.3 

1. 

It 

10/71 

71.6 

17 

II 

II 

II 

61.9 

.8 

II 

73.3 

3. 

II 

III 1/70 

54.7 

94 

.0 

.01 

12.5* 

17.4 

.3 

.0 

59.6 

7.7 

.0 

3/70 

65.6 

53 

1 

II 

II 

50.9 

.7 

II 

70.0 

5.6 

II 

6/70 

80.9 

8 

II 

II 

II 

221. 

2.7 

II 

82.4 

1.3 

II 

10/70 

82.7 

12 

II 

II 

II 

166. 

2.1 

II 

78.7 

2.7 

II 

1/71 

58.8 

92 

II 

II 

II 

19. 

.3 

11 

63.4 

7.6 

II 

3/71 

60.1 

65 

II 

II 

II 

51. 

.8 

II 

65.2 

6.2 

If 

6/71 

71.0 

8 

II 

II 

II 

843. 

10.6 

II 

80.8 

2.5 

II 

10/71 

76.0 

9 

II 

II 

II 

99. 

1.2 

II 

77.5 

2.4 

II 


* null hypothesis can be rejected at a = .05 level 

** I - New Orleans; II - Orlando; III - Daytona Beach 

*** data set consists of dally observation for a monthly interval, only 
these selected months are presented. 


o o o o 
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TABLE 5 

Dally Maximum Wind 


Date*** 

A 

Normal 

A 

A 

Exp 

A 

A 

A 

Gamma 

A 

A 

A 

Extremely 

A 

A 


P 

0 

Q 

X 

Q 

y 

d 

Q 

a 

6 

Q 

1/70 

9.6 

7 

.0 

.10 

11.2* ** 

14.8 

1.5 

.0 

11.8 

2.1 

.0 

3/70 

9.8 

6 

ii 

M 

11 



— 

11.2 

2.0 

ti 

6/70 

6.8 

6 

M 

.14 

9.7* 

7.7 

1.1 

.0 

8.7 

1.8 

ti 

10/70 

7.6 

9 

1.3 

.13 

9.3* 

5.6 

.7 

II 

9.6 

2.3 

it 

1/71 

8.4 

12 

.0 

.11 

8.7* 

4.7 

.5 

II 

10.6 

2.7 

it 

3/71 

9.7 

7 

m 

.10 

11.1* 

11.6 

1.2 

It 

11.6 

2.1 

it 

6/71 

5.3 

2 

i. 

.18 

10.4* 

8.6 

1.6 

II 

6.1 

1.2 

it 

10/71 

4.7 

5 

.3 

.20 

9.1* 

6.5 

1.3 

II 

7.1 

1.6 

ii 

1/70 

9.6 

10 

.0 

.10 

10.4* 

7.8 

.8 

.0 

11.7 

2.4 2 

.4 

3/70 

10.3 

10 

II 

.04 

10.6* 

8.3 

.8 

II 

12.3 

2.4 

.0 

6/70 

8.4 

4 

II 

.12 

11.1* 

14.1 

1.6 

II 

9.8 

1.6 

II 

10/70 

8.8 

6 

(I 

.11 

11.1* 

10.9 

1.2 

It 

10.5 

1.9 

ll 

1/71 

8.P 

7 

II 

.11 

10.7* 

8.7 

.9 

II 

10.4 

2.0 

II 

3/71 

10.1 

11 

M 

.11 

10.7* 

10.9 

1. 

II 

12.7 

2.5 

II 

6/71 

7.4 

3 

It 

.13 

11.3* 

15.6 

2. 

II 

8.5 

1.3 

II 

10/71 

6.8 

5 

II 

.14 

11.0* 

7.1 

1. 

II 

8.2 

1.7 

II 

1/70 

9.2 

5 

.0 

.10 

11.3* 



.... 

10.5 

1.8 

.0 

3/70 

8.8 

6 

II 

It 

11.2* 

11.8 

1.3 

.0 

10.5 

1.9 

II 

6/70 

9.0 

7 

II 

II 

10.9* 

15.6 

1.7 

II 

11.1 

1.9 

II 

10/70 

10.3 

13 

II 

II 

10.3* 

8.6 

.8 

II 

12.8 

2.7 

It 

1/71 

8.0 

7 

II 

II 

II 

8.8 

1.1 

II 

4.4 

2. 

II 

3/71 

9.5 

11 

It 

II 

II 

10.5 

1. 

II 

12.0 

2.4 

II 

6/71 

7.3 

3 

II 

It 

11.5* 

21.9 

2.9 

II 

8.5 

1.2 

II 

10/71 

7.5 

6 

11 

II 

10.7* 

9.7 

1.2 

II 

9.3 

1.8 

II 


* null hypothesis can be rejected at a «.05 level 

** I - New Orleans; II - Orlando; III - Daytona Beach 

*** data set consists of dally observation for a monthly Interval, only 
these selected months are presented. 


'■jlibilii'll* 
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TABLE 6 

Extreme Monthly Temperatures 


Site 

A 

Normal 

A 

Exponential 

A A A 

Gamma 

A 

A 

Extreme 

A 

A 


V 

C 

0 

Q 

X 

Q y 

s 

Q 

a 

B 

Q 

I 

82.5 

1.5 

.0 

.012 

16.9* 


— 

83.3 

.98 

.00 

II 

82.9 

1.3 

.0 

.012 

16.9* 


— 

84.8 

.88 

.00 

III 

81.1 

.8 

.0 

.012 

16.9* 


--- 

81.7 

.67 

.00 


* null hypothesis can be rejected at a * .05 level 

** I - New Orleans; II - Orlando; III - Daytona Beach 


TABLE 7 

Extreme Monthly Winds 


Site 

A 

Normal 

AA 

A 

Exponential 

A A 

A 

Gamma 

A A 

A 

Extreme 

A 

A 


u 

C 

0 

Q 

X 

Q 

V 

6 Q 

a 

6 

Q 

I 

n.i 

17 

.26 

.09 

13.6* 

1.4 

.13 .0 

15.4 

2.9 

.9 

II 

11.2 

11 

.0 

.08 

14.1* 

1.1 

o 

• 

o 

13.6 

2.5 

.6 

III 

10.3 

9 

.0 

.09 

14.5* 

1.6 

.16 .0 

13.0 

2.2 

.4 


* null hypothesis 

** I - New Orleans 

can 
; II 

be rejected at a 
- Orlando; III - 

* .05 level 
Daytona Beach 
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Experiment 3 

In this section the procedure was applied to a data set found In 
Haggard et. al. (1973). In their paper, they analysed a meteoroloolcal 
data set consisting of maximum rainfall amounts In the Appalachian region 
resulting from tropical disturbances. In their paper they satisfactory 
modeled the data set with a Gamma distribution. In this section, I wanted 
to determine If the GMCS procedure would indicate that the Gamma 
distribution would provide a satisfactory fit. Also, since the original 
authors were Interested In making probabilistic forecasts, I have Included 
the slmlllar forecasts based upon the GMCS fitted distribution. The results 
for the test of fit are summar'zed In Table- 7. Table 8 contains a 
comparison of the GMCS fitted Gamma distribution with the results found 
In Haggard et. al. (1964). 


TABLE 7 

GMCS Procedure for Maximum Rainfall within the 
Appalachians 
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Data 

Set** 

A 

Normal 

A 

Exp 

a 

i 

A 

A 

Gamma 

6-1 

A 

Extreme 

/V A 

A 

Haggard et. al . 
Result 

A A 1 


u 

a 2 

Q 

A 

Q 

Y 

Q 

a 

6 

Q 

Y 

6 ” 

A 

7.29 

jJ. 3 

1.75 

.14 

5.14* 

1.9 

3.85 

.14 

16.3 

4.4 

.04 

2.2 

3.33 

B 

8.08 

53.5 

1.24 

.12 

4.70* 

2.2 

3.85 

.09 

16.6 

9.6 

.03 

2.8 

2.88 

C 

9.37 

55.6 

.42 

.10 

3.40* 

1.9 

5.07 

.00 

15.9 

5.2 

.05 

1.9 

4.73 

D 

10.18 

55.3 

.32 

.09 

3.90* 

2.2 

4.56 

.00 

16.6 

5.2 

.03 

2.6 

3.87 


A' 

7.18 

39.7 

1.23 

.13 

5.05* 

2.1 

3.4 

.06 

14.2 

4.0 

.02 

2.2 

3.1 

B' 

7.94 

41.8 

.86 

.12 

4.78* 

2.4 

3.4 

.04 

14.7 

4.2 

.02 

2.9 

2.6 

C* 

9.2 

47.9 

.26 

.10 

3.73* 

1.9 

4.8 

.02 

14.5 

4.9 

.04 

2.0 

4.5 

D' 

10.0 

46.5 

.18 

.09 

4.24* 

2.3 

4.3 

.00 

15.2 

4.9 

.02 

2.7 

VO 

• 

CO 


null hypothesis can be rejected at a = .05 level 

A - maximum 24-hour precipitation all storms. B - maximum 24-hour 
precipitation from no more than one storm per year. C - maximum precipitation 
totals from all storms. D - maximum precipitation totals from no more than 
one storm per year. A' - O' - same as A - 0 except using 27 inches for 
Camille rather than 31 inches. 
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TABLE 8 

Expected Probabilities of Exceeding Arbitrary 
Precipitation Amounts Over the Appalachian Region 


Precipitation Data Sets** 

In Inches 


ABC D 



I* 

II 

I 

II 

I 

II 

I 

II 

1 

.976 

.966 

.992 

.980 

.980 

.980 

.994 

.987 

2 

.909 

.890 

.954 

.926 

.926 

.930 

.968 

.950 

3 

.817 

.797 

.887 

.850 

.850 

.863 

.923 

.894 

4 

.716 

.698 

.801 

.764 

.764 

.788 

.826 

.827 

5 

.615 

.602 

.705 

.674 

.674 

.705 

.792 

.754 

6 

.519 

.513 

.607 

.587 

.587 

.632 

.716 

.680 

7 

.433 

.432 

.513 

.505 

.505 

.559 

.639 

.607 

8 

.357 

.362 

.427 

.431- 

.431 

.490 

.565 

.537 

9 

.292 

.301 

.351 

.364 

.364 

.427 

.494 

.471 

10 

.237 

.248 

.286 

.306 

.306 

.371 

.429 

.412 

15 

.077 

.090 

.091 

.118 

.118 

.172 

.191 

.195 

20 

.023 

.030 

.024 

.042 

.042 

.075 

.077 

.085 

25 

.006 

.009 

.006 

.014 

.014 

.031 

.029 

.036 

30 

.002 

.003 

.001 

.005 

.005 

.013 

.010 

.014 



A' 


B' 

C' 



D' 


I 

II 

I 

II 

I 

II 

I 

II 

1 

.978 

.972 

.993 

.985 

.981 

.977 

.995 

.980 

2 

.913 

.900 

.959 

.934 

.932 

.924 

.971 

.956 

3 

.821 

.806 

.893 

.858 

.865 

.855 

.926 

.903 

4 

.717 

.704 

.806 

.768 

.789 

.779 

.866 

.838 

5 

.613 

.603 

.706 

.673 

.710 

.700 

.794 

.765 

6 

.515 

.510 

.605 

.580 

.631 

.623 

.717 

.690 

7 

.427 

.425 

.507 

.492 

.556 

.500 

.639 

.615 

8 

.349 

.352 

.418 

.413 

.486 

.482 

.562 

.544 

9 

.283 

.288 

.340 

.344 

.423 

.419 

.489 

.477 

10 

.227 

.235 

.273 

.284 

.365 

.364 

.422 

.416 

15 

.070 

.078 

.079 

.098 

.165 

.169 

.182 

.192 

20 

.019 

.023 

.020 

.031 

.070 

.073 

.070 

.082 

25 

.005 

.006 

.002 

.004 

.028 

.031 

.025 

.033 

30 

.001 

.002 

.001 

.003 

.011 

.015 

.008 

.013 


* 1- Haggard et.al. Gamma distribution; II- GMCS Gamma distribution. 

** Same as Table 7 
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Conclusions 

* 

A goodness of fit procedure based upon the theoretical work of 
Dahlya and Gurland [(1970a), (1970b), (1972)] is presented. The proce- 
dure has been documented In the computer software package (Appendix A), 
f Several examples using meteorological data sets are analysed using this 

procedure. The principle advantages of this procedure over existing 
goodness-of-fl t tests lies In the ability to test for several distrl- 
t butlons using a single user setup. This advantage stems from the 

: freedom of testing a distribution without having to specify all the un- 

known parameters of the tabular values of the null distribution. 


References 


Dahiya, R.C. and Gurland, J. (1970a). A test of fit for continuous 

distributions based on generalized minimum chi-square. Statistical 
Papers in Honor of George W. Snedecor , T.A. Bancroft, editor. 

Dahiya, R.C. and Gurland, J. (1970b). Estimating the parameters of a 
gamma distribution. MRC-TR #1067. 

Dahiya, R.C. and Gurland, J. (1972). Goodness of fit tests for the gamma 
and exponential distributions. Technometrics, vol . 14, no. 3, 
pp 791-801. 


107 


Appendix A 

User setup for Gurland's (GMCS) procedure 
JOB CONTROL PARAMETERS 


CARD 

COL 


1 

1-5 

IUNIT 


6-10 

NOB 


15 

ICOR 


20 

IDIST 


25 

30 

35 


2 1-80 NFORMT 


DESCRIPTION 

INPUT DEVICE for DATA. 

Number of observations to be fitted. 
ICOR = 0. 

1 NORMAL distribution fitted. 

0 NORMAL distribution not fitted. 

1 Exponential fitted. 

0 Exponential not fitted. 

1 Gamma distribution fitted. 

0 Gamma distribution not fitted. 

1 Extreme value distribution fitted. 

0 Extreme value distribution not 
fitted. 

Format for input raw data. 

Input raw data 


3+ 
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MAIN 

GCALC 

RHAT1 

RHAT2 

TRIPLE 

AHAT 

QHAT 

GREXTR 

GRNORM 

GREXPO 

GRGAMM 

DGMPRD 

DMIN 


Program Description 

main program; Input job parameters 
calculates the coefficients for matrix G. 

A 

calculates the matrix R for exponential dist. 

A 

calculates the matrix R for other dist. 
calculates matrix product x*y*z. 

a • 

calculates matrix A. 

A 

calculates matrix Q. 

performs goodness of fit for extreme value 
distribution. 

performs goodness of fit for normal dist. 
performs goodness of fit for exponential dist. 
performs goodness of fit for gamma dist. 

IBM matrix multiplication 
IBM matrix inversion 


Subroutines Needed By A Given Routine 


MAIN - GRNORM, GREXPO, GREXTR, GRGAMM 

GCALC 

RHATi - DGMPRD 

RHAT2 - DGMPRD, DMINV 

TRIPLE - DGMPRD 

AHAT - DGMPRD 

QHAT - DGMPRD 

GREXTR - GCALC, TRIPLE, DMINV, RHAT 1, AHAT, QHAT, DGMPRD 

GRNORM - GCALC, TRIPLE, DMINV, RHAT 2, AHAT, QHAT, DGMPRD 

GREXPO - Same as GREXTR 


GRGAMM 


Same as GRNORM 


FURTRAN IV G LEVEL 21 


MAIM 


DATE - 78 1 92 



0001 

0002 


0003 

0004 
0003 
000b 
0007 
0006 

0009 

0010 

0011 

0012 

0013 

001 4 

0015 
OOlo 

0017 

0018 
GO 19 

02 0 ' 
021 
022 
02 3 

024 

025 


b 

0027 

0028 
0029 
003 0 
0031 


32 

33 

34 
003 5 
0036' 
0037 


0038 

003 9 

0040 

0041 
0042' 

0043 

0044 

0045 

004 6 

0047 

0048 

0049 

005 0 
0051 
00 52 

0053 

0054 
0065 


9030 

1 


IMPL I C IT RE AL*8 ( A-H * Q-Z) 

DIMENSION RAW (81 * CUNL ( 8 ) * CENRL ( 8 ) , G( 4* 4 ) * X( 1000)* ID 1ST (10), 
i' NFURMT(20)»XJll4#4),EX0B( iOOOl 

DIMENSION L I ML ( 33 ) 

COMMON / MOMENT / R A W* CUML* CfcNRL* S» NOB '* 

COMMON / NUMBER/ XO I V > X ME AM# X VAR , X GE OM, l JM I T* I COR# P I » STD 
READ(5*1*£N0«9999) I UNIT# NOB* ( I 01 ST ( 1 i * I >1* 5 ) 

FORMAT (915) 

RE A0( 5*2) (NFJRMT( I ) »I»1*20V 

FORMAT ( 20A4 ) 

IFU01SM1I .EG. 5) GO TO 9004 
READ! IUN1T*NFURMT) (X(J)» J > 1*N0B) 

DO 152 I»1*N0B*12 

XMAX « X m 

xm in « xi i i 

K » I 4l 

L « 14-11 

DO 151 J > K*L 

IF (XI J) .GT. XMAXI X M A X » X ( J I 

IF (X(J) . L T . XMIN) XM IN * X ( J ) 

CONTIMJE 

WRIT£(6. 153) XMAX* XMIN 
FORM AT ( T2 5* 5F5.1 ) 

ICHECK « 0 

DO ill J * 1»NQB 

IF ( X( J I .LE. 0.) ICHECK « 1 


WRIT£(6» 125) 

WRITE! b* 123) ( X( J ) » J * 1* NOB ) 


9031 


XMIN 


ICHECK 


XOIV * 
XMEAN - 
X V A < « 

SVAR * 


S 
X 

XM4 - 0, 
SM2j^ 0, 
SM3 * 0, 
SM4 « 0, 


0FL3AT ( NOB ) 

0.0 


a 


Quality 


9032 

c‘ 

c 


PI « 3.1415926 
SOIV « XOIV - 1. 

DO" 9 CO 1 T ■ ” 1* NOB 

XMEAN • XMEAM ♦ X(I) / XOIV 

SUM « S JM * DABS ( X ( 1 ) ) 

IF (SUM .LE. 0. ) SUM»0. 1 

SO « OLOG(SUM) / XDIV 

XGE3M - DEXP(SD) 

CONTINUE " 

00 9 CO 2 I « 1* NOB 

XVAR « X VAR *■ ( X(I) - XMEAN )**2 / XDIV ~ 

XM3 _ * XM3_ ♦ ( X(II - XMEAM ) _** 3 / XDJV^ 

XM 4 - XM4 «• (XII) - XMEAN I ** 4 / XDIV 

SM2 ■ SM2 ♦ X ( I ) **2 / XOIV 

SM3 » SM3 ♦ X ( I ) **3 / XOIV 

SM4 « SM4 4 X( I )**4 / XDIV 

STD » DSURT(XVAR) 

CONTINUE 




c 

c 

c 


00 2b 
0027 
0020 
0020 


OObO 

OObl 


0 Obi 
00b3 


0 Ob* 
0065 
OObb 
0067 


0060 

0069 


LOOP TILL ALL DISTRIBUTION REQUESTS HAVE BEEN SATISFIED 


DO 9003 1 - :.a 

IF umST(I) .IE. 01 .OR. UDISTCU .GT. A ) ) GO TO 9003 
I0U1 • IDtSTC I ) 

GO TO ( 11. 1 1.13.1A). IDUM 


C 

C 

C 

11 

C 

C 

C ' 

12 

c _ 
c 

_c_ 

13' 


NORMAL 

CALL GRN CRM( XM3 ) 
GO TO 9003 

EXPONENT! AL 


ID I ST * L 


1 1) I ST 


CALL GREXPO! SM2.SM3.SMA.X) 
GO TO 9003 


GAMMA 


121 

C 

j: 

C 

LA 

9003 


IF C ICHECK . EO . 0 ) 

WRITE (2.121 ) ICHECK 
FORMAT! 10X. 251 I2.1X)) 
GO TO 9003 

EXTREME VALUE 


IOIST « 3 

c a ll grgavTFx 7 sTT 2 .1 mb7s~ma F 


ID I ST « A 


CALL GR5XTR|XI_ 
CONTINUE 



C 

C 

BIVARIATE NORMAL IDISTI1) * 5 

0071 

C 

900 A 

RE AD ( 1JNIT.NFQRMT ) I ( X ( ( J-l ) *2 + 1 ). X( ( J - 1 ) *2 + 2 ) ) . J-l.NOBI 

00 72 
0073 

123 

CALL 81VARIK.NUB. IUNIT) 
FORMAT ( T 25. 5F12.5) 

00 7 A 

122 

C 

FORMAT! 1 F1///1H0.T51. * THE OBSERVATIONS*.//) 

0075 

007b 

9999 

25 

MRlTElb. 25 ) 
FORMAT!* 1‘ ) 

0077 

0078 

10 

REWIND 9 

REAIM9. 15. END-20) LINE 

0079 

0080 

15 

FORM AT ! 3 3AA ) 
WRITElb.15) LINE 

0081 

0082 

0083 

20 

GU TO 10 
STU> 

END 
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FORTRAN IV G LEVEL 21 OCALC OATE ■ 78192 

SUBROUTINE SCALCIICOR)" 

CALCJLAT 6'T‘FURTfRSf FOUR 01 SH I BUTTONS 
IMPLICIT REAL*8 IA-H * Q-2 ) 

DIMENSION R4WI 0) » Gt 4* 4) *CUML IBI*CENRH8)*AI1000)*B(1000) 
COMMON /MOMENT/ RAW»CUML*CfcNRL* 3* NUB 

XN ■ OFLUATtMOB) . . _ . 

00 100 I - 1*4 
00 ICO J > 1*4 

Gl 1 * J ) « R AR I I ♦ J ) - RAW! I )*RAWl Jl 
CUNTINJE ' 

RETJRN 

END _ 


FORTRAN IV G LEVEL 2i RHATl DATE « 78192 

0001 


0002 

0003 

0004 

0005 

0006 

0007 

0008 

0009 

0010 
0011 


SUBROUTINE RHAT 1 1 W* S IG I * R ) 

C 

C CALCJLAT E VECTOR R HAT FOR EXPONENTIAL DISTRIBUTION 

C 

IMPLICIT REAL*8 CA-H * Q-Z) 

U1MENSIJN WI4)*$IGII4»4)»RI4*4)»0UMI4)*XI 1)*FJURI4»4) 
CALL OGMPROI •»* S IG I * JUM * l » 4 * 4 ) 

CALL OGMPROI OUM»W» X* 1*4* l ) 

xu i « i.o 7 xt D 

CALL DGNPRDI X*rt*DJM» l# 1*4) 

CALL OGMPROI W* DUM * F0UR*4 * 1 * 4 ) 

CALL OGMPROI FOUR* S I G I * R* 4* 4* 4 ) 

RETJRN 

END 


0001 

C 

* C 

c 

0002 

0003 

0004 

0005 

0006 

0007 

0008 

0009 1 03 

0010 
0011 


FORTRAN TV "6 LEVEL 21 RHAT2 DATE « 78192 I 

0001 ” SUB ROUT I NE RHAT2 1 W# S IG I * R ) 

R HA r MAT FfxYTx 2T ~F~U R GAMMA* NEGTl N* NORMAL” 

0002 ~ IMPL K IT REAL** Ya'-H 7 0-Z') 

000 3 DIMES SU N W 14* 2) * S I G 1 14*4 ) * R ( 4 * 4 ) * W TI2 _» 4 ) *JKJM ( 2* 4)*X|2 * 2 ) » 

£ F UU R 1 4 * 4 ) * M 1 2 ) * L 1 2 ) 

0004 00 900 0 1 « i*2 

0005 

0006 _ 9000 

0007 

0008 

0009 " 

0010 

“00 U 

0012 

0013 ’ 

0014 


DO 9000 J ■ 1*4 

WTI 1 * J ) - WI J*I) 

CALL DGMPKCI WT* SIGI *DUM#2*4*4) 

CALL DGMPRDt OUM»W» X* 2* 4*2) 

CALL” DM I NV (X*2*DET*L»M) 

CALL DGMPRDIX*WT»DUM* 2*2*4) 

C A L L* D GM PRC l W » D UM * F J U R * 4 * 2 *7 ) 
CALL OGMPROI FUU F* S I G I * R* 4 * 4* 4 ) 
RETJRN 
END 


JRTRAN IV 6 LEVEL 21 


DATE 


113 



A HAT 


0001 


0002 

0003 

000A 

0005 

000b 

0007 

oooa 

} 0009 

0010 
0011 
0012 


“SUBROUTINE AHATI S I Gl * R# A I 

c 

C CALCULAT i A HAT 

C 

IMPLICIT REAL** < A-H * 0-21 

DIMENSION SIGItA#M»pRtA#AI#AtA*A»»_Rl(A**l 

ooT i ■ iVa 

DO 1 J • 1 * A 

RI(f#J) ■ -RU#J) 

IF ( 1 .NE. J I GO TO 1 
Kl(l»J) - RIU»JI ♦ 1.0 
1 CONTINUE 

CALL DGNPROI S IGI » R l » A» A* A# A ) 

RETURN 

END' 


: ORTRaN IV G LEVEL 21 TRIPLE DATE • 7*192 


i 

> 0001 


0002 ' 
0003 
000* 
0005 
00 ’ 06 

0007 

0008 

0009 

0010 


SUBROUTINE T RIPLE ( X. Y» Z I 
C 

C CALCULATE X * Y * X TRANSPOSED AND RETURN VALUE IN Z 
C _ _ 

IMPLiC IT REAL*8 U-H » 0-Z) 

DIMENSION X(A*Al#Y(A»AI»Z(A>AI.UUMlA,A). XTI A» A I _ 
DO 1 I « 1*A 

DO 1 J « 1*A _ 

1 XT ( I # J 1 * X( J » 1 ) 

CALL DGNPRQ(X#Y»DUM»A»A»A) _ 

CALL DGN PR C( DUM»XT#Z»A*A»A) 

RETURN 

END 


: ORTRAN IV G LEVEL 21 OH AT 


0001 


WO 2 
0003 
000 A 
JD005_ 
0006 

0007 

0008 


SUBROUTI NE QHAT ( XN» H. A.Q ) 

C 

C CALCULAT L CHI-SQUARE Q HAT' 

C 

"IMPLICIT REAL** I" A-H / 0-Z ) 

DIMENSION H ( A I # A ( 4 » A I # D U M ( A ) j X X (_1 1_ 

CALL DGN PRDTh/a# DUM» 1 # A. A I 

CALL DGNPRD(DUM#H»XX*L#4# II 
Q * X X (LI * X N * 

RETURN 

END "" ' 


Original 

of , 


page is 


i 
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JRTRaM IV 
jOOL 


-602 

-*003 


004 

-005 

-006 

:007 

.008 


-009 

-JlO 

-Oil 

.012 

-013 

014 

- Oi 5 
-016 

017 

-018 

- Ji 9 
-020 
-021 
-022 
-023 
-024 
-02 5 


-026 

-027 

-028 

0029 

-030 

-031 

-032 

-033 
- 03 * 
-035 
- 0 3 fa 

-037 


-038 


C LE/tL 21 GREXTR OATE ■ 78192 14/47, 

SUBROUTINE GRfcXTRtXI 

C_ . 

C GURLAMO ROUTINE FOR EXTREME VALJe JlsTRlBUTIUN 
C 

IMPLICIT ReAL*8 tA-H#0-Z) 

01 HEMS ION XJl(<t»4l»RAWm»CUm.U)»CENRLm»G(4»4)»tfl4i#HU)» 

£ JUMt 4#4l»SlGl(4*4l#Lt4l*Mt4)#THETAlA)*Rt4*4l#Al4#4|# 

C _ XtlCOOHBriATtiOG)»TOENOMtlOO) 

C J M M ON / MO M S N T / R A W» CUML # C ENRL * G» NOS 

CUMMJN / NUMBER/ XO IV#X PEAN*XVAR» XGEOM# IUMIT» ICOR»P 1>ST0 
C 

XN • DFLOATtMUB) 

ZERO* 0. 0 

ONE- 1.0 _ 

' C ’ 

C CALC JL AT E EXTRtME CUHULANT MOMEMTS 

C 

C POT COMLU-ll IN PLACE OF COHLlII 14 ORDER TJ MAXc THE SAME 

C SUBSCRIPTS OF H-VECTOR AS THAT OF THE JACU8IAN MATRIX 

_C _ _ . . 

eVum - o7o”~ 

SIX-b. 

BETA • OSQRT ( SI XI * STO / PI • 

8 • BETA 

DO 1 I • 1 » NOB 

1 . esum_ - j sum *j)Exp( x<Dj/bi_ „ 

' AL p HA • B * OL 0 Gt E S J M f - B * DLOGtXDIV) 

EMEAN - ALPHA - 0.577216*8 
EM03E ■ ALPHA 

EVAR ■ * l ** 2 * 8 ** 2 / 6. 

CuNLtl) • 1.6*5*B**2. 

CUML ( 2 ) - 2 . 396*8 ** 3_. 

CUML 13) - '6'.494*B**4. 

CUMLt 4 ) « 24.886*B**5. 

CUML ( 5 I • 122.078*8**6. 

CUML (6) - 72 6. 01*B**7 . 

CUMLt 7) « 5060.5*15*8**8. 

__C 

Ci ■ CUML ( 1 ) 

C2 • CUM 1121 

C3 CUMLt 31 "" ' 

C4 ■ CUMLt 4 1 

C5 • CUMLt 51 

C6 - CUMLt 6 1 

C7 ■ CUMLt 7 1 

C _ _ 

RAWt 1 1 - XME AN 

RAW I 2) • Cl ♦ XME AN** 2 __ _ 

RAW (3 ) ■ C2 * 3 . *C 1*XMEAM ♦ XMEAM**3 
RAW 1 4 1 ■ C 3 ♦ H.*C2*XMEAN ♦ 3.*C1**2 * 6. *C_1*XMEAN**2 
t ~ * XME AN* *4 

RAWt 5 ) « CA ♦ 5.*C3*XMEAN ♦ 10.*C2*C1 ♦ 10. *C2*XMEAN**2 
t * 15,*C1**2 * XME AN * 10. *C i*XME AN* *3 ♦ XM£AN**5 

C 

RAWt 61 ■ C 5 * 6, *C4*XMEAN ♦ 15.*C3*C1 ♦ 1 5. *C 3*XME AM** 2 ‘ ♦ 

t 10. *C2**2 * 60.*C2*Cl*XME AM ♦ 20. *C2*XMEAN** 3 

t *"li.*Cl**3 ♦ 45.*Ci**2 *XMEAM**2 ♦ 1 5. *C l * X MEAN"** A * 

L XMEAN**6 


-MGWAL PAGt li 
* POOR QUALITY 




i 


JRTRAR IV C LE/tL 21 
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CREXTR 


DATE ■ 74142 


14/ 47/ 


0039 


RAMt 7 ) 


£ 

e 

£ 

t 


J040 


RAMI 8 1 


£ 

& 

£ 

& 

£ 

£ 

£ 


Co ♦ 7.4C5*XMEAN ♦ 21.« , C4*:i ♦ 21.*C4 *XMEan * *c 

♦ 35\*C3*C2 ♦ i06.*C3*C.4XMEAN ♦ 35 . *C 3*XMEAN**3 

♦ 70.*C2**2 *XMEAN ♦ 135.*C2*C1**2 ♦ 210.*; 2 

*Cl*XM£AN**2 ♦ 35.*C2*XMEA1**4 ♦ 105,*Ci**3 *XMEAN 

♦ 105. *11**2 *XMEAN**3 ♦ 2&«*C1*XRcAn** 5 * XMfcAU**? 

C7 ♦ 8.*C6*XM£AN ♦ 28. *C5*Ci ♦_ 24 . *t 5 *XMfc AM**2 ♦ 

W.*'C4*C2 ♦ 1 b8 . *C ^**C 1 ♦ <fit A 4 » 36.*C**XMEAN**J ♦ 

35.*C3**2 ♦ 280.*C3*C2*XMeAN ♦ 210. *C 3*C a**2 ♦ 

423 .*C3*Cl*XMEAN**2 ♦ 70. *C3*X1EAN**4 ♦ 25Q.*C2**2 *C1 

283.*C2**2 *XMEAM**2 ♦ 443. *C2*C 1**2 *XMcAI ♦ 5 o,j.*C 2*C 
*XMEAN**3 ♦ 56.*C2*XMcAN**5 ♦ 1 05 . *C 1 ^ % 

♦ 420.*Li**3 *XMEAN**2 ♦ 2i0.*:i**2 *XM£AN**v 

♦ ,28.*Cl*XMfcAN**& *• Xit AN** 8 


.041 

r 

CALL 

GCAL Cl 1C0RI 


w 

c 

r 

INITIALIZE w 

>042 

w 

Ml 11 

■ • 4 J6 

.04 3 


4(1) 

« ^ .713 

.044 


Ml A 1 

• 3.853 

04 5 


Ml 4) 

•4. 906 


C 

c 

IMiriALIZE H 

.046 


Hi 1 1 

• C UML I 2 ) / CJML 1 1 ) 

0.4 7 


HI 2 1 

• Z UML I 3 J / CUML I 2 ) 

.048 


Ht 0) 

« : UML C 4) /CUML 1 3) 

<049 


HI 4 ) 

> ; UMLI 51 /CUMLI 4) 


C 

c 

c 


INITIAL 1 Z £ 


J1 


-350 

.051 

352 

.053 

>054 

.355 

.056 

.057 

.058 

.059 

.060 

.061 


062 
063 
0t>4 
. 365 
.366 
.067 


TO 120 


12) 


C 

C 

C 


00 120 2 >1/4 
00 120 J *i » 4 

iFici.to.j i.OR.m-n.Eo.jii gj 

XJ1( 1. J) *ZfcRO 

CONTINUE 

XJ 111,1} - ONE 

XJ1I2,2) - l./CUMLlil 

X i 1 1 2 , 1 ) « -CUMLC 2 l/CUML I 11**2. 

XJ 1 1 3, 3 ' • 1./CUNLI2I 

* -CUMLt 4 J/CUMLJ 3J*_?2_.„ 
Xjlii.2) • -CUMLI3J/CUMLI2I**2. 
XJ 1 ( 4 , 4 ) » l./CUMLlil 


CALCULATE CHI-SQJARE TEST ANO EXTREME PARAMETER 


CALL TRIPLEIXJ1,G,SIG1 I 
CALL DMINV (SIGI#<.»0ET»L»M) 
CALL RHA T1(M,S1G1,K) 

CALL AHA T(SIG1,R,A} 

CALL OHA Tl XN, H, A, U ) 

CALL DG1PR0(R,H, THETA, 4, 4. 1) 


. Ob 8 
.0 69 


iRfrrrb. i2 5f 

MKlTE(6, 1221 


IXIJI, J-1#N0B) 
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JRTRan' lV~G Lt/eL 21 " GREXTR OATE » 70192 l%/*7/ 


.370 

WRITE!©. 123) 

emean 

j71 

WRITEI6* 12%) 

STD 

372 

WRITE To* 12b) 

EMJDE 

.373 

WRITE 16* 127) 

EVAR 

37% 

WRITE! o»i2"3) 

XV A R 

375 

WRITElo* 12 9) 

XMEAN 


37b" WRIT£f9* 130)' AL PHA* 8*9 

377 WRITE ! 5*121) ALPHA*8*Q 

IT* TTi FORMAT!// ,T25* • PARAMETERS : ALPHA* • * E 15 . 5 * 10 X* ' geTA* '*£15.5* 

£ //*T39*' ♦ ♦♦! CHI-SQUARE VALUE )*** **£15.5) 

379 * 122 FORMAT ( X J5* SFTO.'S)"' 


opo 123 F0RMAT!///*T37*'THE MEAN OF THE EXTRtME VALUES IS**Fi5.7#/| 

381 127 FORMAT)/ //* T37* 'VARIANCE OF THE EXTREME VALUES ISSF15.7,/) 

382 12% F0R1ATC/// »T37 **THE STANOARO DEVIATION I S * # 8 X* F 1 5. 7# / I 

‘j 83 125 "FORMAT ( 1 Hi / / /TrtO* T 5 i * ' f HE OBSERVATi J nS' 7 lH* 

£ T%1 * * GJR LAN3S PROCEDURE FOR EXTREME VALUES'*//) 

38% 12b F0RM/T(///*r37*'THE MODE OF THE EXTREME VALUES I S' * F 1 0 . 7* / ) 

385 125 F0RMAT(///*T37*'THE SAMPLE VARIANCE I S' » 1 IX* F 15. 7) 

C INITIAL IZE M 

j 86 i_2 9 FORMAT !/// * T37*'THE SAMPLE MEAN 1 S ' » 1 5 X« F 15. 7 ) 

387 133 FORMAT(//*T25*' EXTREME' PARAMETERS: AL P H A » • # F b . 2* 5 X * • BETA*'** 

£ F 6.2* //* T39* * ♦ ♦»! Crll-SQUARE VALUE )•** '»F13.3i _ _ 

388 RET JRM ‘ 

389 END 


: JRTRaN 

0001 


0002 

0003 


000 + 
JO 05 
000a 

0007 

0008 
3009 
0010 


oo n 
0012 
0013 
aO 1 + 
0015 
JO i 6 
00*7 
j0*8 
.019 


0020 

j'021 


.022 
.023 
j02 ** 
0025 


0026 

.02 7~ 

0028 

.029 

j030 

0031 

.032 

-03 3 

.03 + 

.036 

.036 

.037 


IV 0 LEVEL 21 
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GREXPQ 


DATE • 78192 1+/+ 


SUBROUTINE 5 REXPUISN 2 . SM 3 »SM+»X) 

C 

C GURL ANl) KJUT INE FUR E\PQNENTIAL DISTRIBUTION 
C 

INPLICIT RfcAL '6 ( A-N * Q-Z) 

JlflENSI JN X 31 ( + * + ) #RAWl 8 ) #CUHL I 8 1 * CfcNRL I 8 )# Gt + * + I # WI ■» I » 
e •h + i»oum( + ' + ) #S 1 G 1 I +# + )#li •*) » ni + )» theta i + i» 

& 4 (+»+:»Ai+»+>»x(iooo) 

COfl'IUI 7 NUN 3 EK / XO I V# X PE AN# XVAK* XGE ON# IUNIT# ICOR#Pl#STO 
CUniJr< / MOMENT / RAW#CUPL#C 6 NRL# 5 #NU 8 

wri tei o» loco 1 1 un i r 

*ooo furiati///#* exponential oistrisution with data from jnit 

& 13 #//) 

XN • OPLLATINOSI _ 

ZERO • 0.0 
ONE • 1.0 
C 

c calcjlate exponential numents 

c 

_ RAWl 1 1 • XMi AN 

RAW I 2 ) ■ 2 * X H E A N * • 2 
RAWI 3 ) ■ 6 * XHEAN **3 
RAWI +1 ■ 2 + * XHfc AN**+ 

RAWI 5 I • 120 * XMEANP *5 
RAWlo) ■ 720 * XMEAN *+6 

RAWI 7 I ■ 5 C + 0 * XHE Af 4**7 _ _ 

RAW! 81 i +C 320 * XMEAN^S 
CALL GCALCUCORI 
C 

C IN 1 TI AL HE W 

c 

00 9000 1 - 1 # + 

9000 till • 1 
C 

c initialize h 

c 

Hill ■ RAW(I) 

Raw’ll * 

HI 31 ■ S H 3 * / SM 2 
HI + ) • SM+ / sh: 
c 

C INITIALIZE J 1 

C 

DO 90 0 1 1 « 1 # + 

00 9001 J • 1 # + 

IF 1 II .EQ. J> .OR. I ( 1 - 1 ) .EO. JIJ GO TO 9001 
XJ 1 I 1 #J) - ZERO 

9031 CONTINUE __ 


' XJ1IMI 
XJ1I2#2) 

ONE 

1.0 / RAWI 1 ) 



XJ1I 3*3) 

1.0 / A AW I 2 ) 



XJl(+#+) 

1.0 / RAW I 3 ) 



XJ i I 2. * ) 

-RAW 12) / RAW(1>**2 



XJ1I 3.2) 

-RAW ( 3 1 / RAWI2)**2 



XJll+#3) 

-RAWI+) / RAW! 3 )**2 




C _ 

C CALCJLATE CHI-SQUARE TEST AND EXPONENT I AITpaRaHETER 

C 


* *** S2& 






mmmm 
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ORTRAM IV G LEVEL El 


GRGAMM DATE - 78192 14/4 


0001 SOtJROUriNt GRGAMMIX»5M2*SM3»SM4) 

C __ 

£ GURL ANO KuUr INE FOR GAMMA DISTRIBUTION 

C 

0002 ‘ IMPLICIT REAL*8 IA-H » J-Z) 

OOoi DIMENSION X Jil4,4)»XJ2l4»4),RAMId),CUMLl8)»C2NRLla)* 

£ 3(4,4i,M(4»2),H(4),DJM(4,4),SlGIt4,4),Lt 4)»Ml4)» 

6 THETA! 4)»RI4,4),A( 4# A ) »X( 1000) » SCO ML (1001 

000 A ' COMMON /MOMENT/ K AW » CUNL # C fcNRL, G »NJ 3 

0005 COMMON / NUMBER / XDIV*XMEAN,XVAR*XGEOM,lUNlT*lCQR#PI,$TD 

0006 WRITEtb, 10CO HUNIT 

0007 1030 FORMAT!///# 1 GAMMA DISTRIBUTION WITH DATA FROM UNIT *,I3»//) 

0008 XN « OFL LiAT ( NUB I 

000 9 ZE RO » 3 .0 

00i,0 UNc * i. C 

0011 IF I XMEAN .LE. 0.1 wR I T£ I b» 21 ) 

0012 21 FORMAT ( * ***NcGATIVE VALUES WRUNG 0 l S TR I BUT I ON*** •» 


C 

C CALCULATE GAMMA MOMENTS 

C 

00x3 ' AX « 0 A 3 S I XMEAN) 

0014 Y v « DL0G10IAX / XGEOM>_ 

0015 YYY « OABS IYY) 

00X6 T 1 « .2500 * 1 1.0 / Y Y Y I * II. 0 «■ DSURTtl. ♦ 1.3333333333 00 *YYY I i 

0017 ' T2 « AX / Tl 

0018 CUMLIi) * T l * T 2 ___ 

0019 ‘ CUML (2 1 « Tl * T 2**2 

OOtO CUML I 3 ) * Tl ♦ T2** , 3 ♦ 2.0 _ 

Q021 CUML l 4 ) * Tl * T2**4 * 6.0 

0022 DO 351 1-1,8 _ 

0023 TL I M IT * 2. *♦(-2511 

0024 XX » OFL OAT I I ) 

0025 IF lITi.GE. 57.57) .OR, C fi" . L E . TL I M I T )) GO TO 98 

0026 RAW! I) « DGAMMAITI *■ XX) / DGAMMA(Tl) * T2**XX 

0027 GO TO 351 

0028 93 RAW! I) * ITUXX-i.) * T2 ** XX 

0029 351 CONTINUE 

C 

C " CALCULATE SAMPLE CUMULANTS FOR GAMMA " 

C 

0030 SCUML(l) - RAWU) 

0031 SCUMLI2) * S.M2 - RAWll>**2 

0032 SCUMLI3) ■ SM3 - 3. *SM2*RAWI 1 ) ♦ 2,*RAW(l)**3 

0033 SCUMLI4) - SM4 - 4 . »SM3»R AW 1 1 ) - 3.*SM2**2 ♦ 12.*SM2*RAWI 1 ) **2 

'fc ' -6 . *RAW 1 1 ) 4*4 ' 

C _ 

0034 CALL GCA LC 1 1 COR) 

C __ 

c 

0035 Wti»i) • ONE 

003b MI 2, 1) « Z ERO ' “ 

0037 W(3»l) • ZERO 

0038 MI 4 » 1 ) * ZERO 

0039 MI 1 » 2 ) * ZERO 

0040 Ml 2, 2 ) « ONE 

0041 Ml 3. 2 ) ■ 2.0 

0042 MI 4,2) « 3.0 * 

C 
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FUft'TAAH IV G LEVEL GRGANN DATE - 78192 14/41 


C " INITI AL I IE H 

C 

JO A3" rtiir- SCUHLtl) 

004* HI2I • SCUMLI2I / SCUNLILI 

00*5 H( 3) • SCUMLI3I / SCJMH 2) 

0046 HI 4 1 • SCUNLI*) / SCUHLIJJ 

"• ‘ c 

C INITIALIZE Jl 

c 

0047 00 LOO I - l# 4 

0048 DO 100 J • 1.4 

0049 IE I . .GT , J I GO TO LOO 

0050 IF (I . E Q. J) XJ1I1.JI - ONE 

0051 IF II .N S. J ) XJit I. J) ■ ZERO __ __ 

0052 iOO CONTI NOE 

0053 XJll2.il - -2 * RAWtll 

0054 XJH3.il * - 3*R Ad I 2 1 ♦ 6*RArtll) 

0055 XJ1I4.1I « -4*K Art I 3 1 ♦ 1 2*RArt l 3 I *RArt 111 - 24*RArt 1 1 1 ♦* 3 

0056 X J 1 1 3 . 2 1 « - 3*R AW I 1 1 

0057 X J1I4.2I « -6*R AW( 2 I » 12»RArtll)»*2 

0058 XJ1I4.3I « -4*RArtlll 

C ______ 

C * INITIAL I ZE J2 

_c __ 

005 9 DO 1 01 T » 1. 4 

0060 DO 101 J * 1.4 

0061 IF I II .fcO. Jl .OR. III-ir.EO. Jl 1 GO TO 101 

0&62 X J 2 f I. J) - ZERO 

0063 ' 101 CONTINUE 

0Q64 X J 2( 1.1) * DNE 

0065 " X J 2 1 2 . 2 1 » 1.0 / CUNLIL) 

0066 X J 2 f 3.3) * L.O / CUNL ( 2 I 

0067 XJ2U.4I - 1.0 / CUNL 13) 

0068 XJ2I2.1I » -CUMLI2) / CUNL I 11**2 

0069 XJ213.2) * -CUML(3) / CUNU2)**2 

0070 XJ2I4.3) » -CUMLI4) / CUMLI3)**2 

C 

C CAICJIATE CHI-SOJARE TEST AND GANNA PARANETERS 

c ' 


0071 

0072 

0073 

0074 

0075 

0076 

0077 
00*0 

0079 

0080 
0081 
0082 

0083 


CALL TRI PLEIXJL.G. OJN) 

CALL TRI PLEIXJ2.DUM.SIGI I 
CALL DMINVISIG1.4.0ET.L.MI 
CALL RHAT21W.SI CI.RI 
CALL AHA T( SIGI.R. A) 

CALL OH A T I X N »H, A. 0 I 

CALL DGNPRDIR.H. THETA. 4. 4. 11 

XR « THcTAli I / THETAI2) 

XL » L.O / THET A( 2) 

MR I T El b. 12 31 XK.XL.Q 

WRITEI9. 124) XR.XL.Q _ 

123 FORMAT!//. T25.‘ PARAMETERS s" k '»' • »E l 5. 5. 10X. • L ANOA >• . E 1 5 . 51 

£ / /. T 39. ' *♦*( CHI-SQUARE VALUE >♦** • » E 1 5 . 5 ) 

124 FORMAT!/ /. T25. • GAMMA PARAMETERS: R- *.F6.2.5X. # L AMD A ■ *.Fb.2» 

£ / /. T39# ' ♦ *♦! CHI-SQUARE VALUE I*** *»F10.3) 


0084 RETURN 

0085 END 
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-URTRaN IV G LEVEL 21 


GRNQRM OATE • 78192 14/47 < 


UOOl SUBROUTINE 5RNGRM IXM3#XI j 

C GURLA^T^QRm 

c . ..... i 

0002 IMPLICIT R£4L*8 (A-H * O-ZI 

0003 DIMENSION X J 1 ( 4* 4 1 » X J2( 4#4 ) » RAW I 8) #C UNL I 8 ) * CENRL 1 8 »• 1 

& 3(4»4)»W(4>2)»H(4)*00M<4'4)»SiGIt4»4)»L(4)»N(4)' 

& THETA(4)#R(4»4)»A(4»4)#X(1 000 1 _ j 

0004 £OMiJN /NUMBER/ XO I V, X PE AN# XVAR - # XGEOM# I UNI T* ICOR* P 1 # S TO 

0005 COMMON /MOMENT/ RAW#CUML*C£NRL* 3* NOS __ | 

000b WRlTE(b» 10C0) 1UNIT " 1 

0007 1330 FORMAT! / //# • NORMAL DISTRIBUTION WITH DATA FROM UNIT *#13*//I_J 

0008 XN ■ OFL OAT! NUB ) j 

0009 ZERO « 3.0 

0010 ONE » 1. C 

C _ _■ 

C CALCJLATE NORMAL MOMENTS 

C _ 5 

0011 CENRH1I • ZERO 

0012 CENRL ( 2) ■ X VAR | 

0013 CENRL i 3 1 » ZERO 

00 1 4 CENRL (41 « 3 * X VAR *42 

00 i 5 ' RAM 1 1 ) * XMEAN 

001b RAW! 21 ■ XVAR ♦ XMEAN**2 

30 £ 7 RAW 13) « 3 * XMEAN * >VAR ♦ X Me AN** 3 

0018 RAW ! 4 ) ■ 3 * XVAR**2 «• b * XME AN**2 * XVAR ♦_ XMEAN** 4 

"0cTl9 RAW! >) « 15 *“XVAR**2 * XMEAN ♦ 10 * XVAR * XMEAN**3 ♦ XME AN**5 

0920 RAW ( b ) * 1 5*XV AR**3 ♦ 45*X V AR* XM t AN ** 2 ♦ 1 5* XV AR*XME AN* *4 ♦ 

£ XME AN**b 

00*1 RAW! 7 ) • 105*XMEAN*XVAR**3 ♦ 84* XV A R* *2 * XME AN** 3 ♦ 

£ 21*XVAR*XMtAN**5 XM£AN**7 

0022 RAW! 8 ) « 105*XV Ak***» ♦_ 4 20* XV AR* * 3* XME AN * *2 f 

£ 2 1 3*X V AR**2*X ME AN** 4 ♦" 28*XVAR *XME AN**b ♦ XMtAN**8 


0023 CALL GCALC(ICQR) 

C 

C IN I T I AL IZE W 

C 


0024 

W( 1* 1) 

* 

ONE 

0025 

W( 2, 1 > 

■ 

ZERO 

002b 

W ( 3 » 1 ) 

S 

ZERO 

0027 

W( 4 » 1 ) 

■ 

ZERO 

0028 

W( 1»2I 

* 

ZERO 

0029 

W( 2*2) 

* 

ONE 

0030 

W!3»2) 

m 

ZERO 

0031 

W ( 4 » 2 ) 

m 

2.3 



C 

INITIALIZE m 


C 


0032 


HID « XMEAN 

0033 


H( 2 I ■ ) LOG 1 CENRL ( 2 ) ) 

0034 


H(i) « X M3 

0035 

C 

H( 4 1 > 3 LOG (CENRL ( 4 ) / 3.0) 


C 

INITI AL IZE Jl 


C 


003b 


00 ICO I - U4 

0037 


DO 1 GO J - l » 4 


0038 


IF I 1 • 3 T. j I GO TO 100 
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-uRTRAM 

U03 9 
J040_ 
0041 
0*042 

0043 

0044 

0045 
004b 
0047 


0048 

0049 
j050 

0051 

0052 

0053 

0054 
00 55 
005 b 

0057 

0058 


0059 
CObO 
00b 1 
00b2 
00 b3 
O0b4 
00 b 5 
00b b 
00b7 
J068 
0069 
_0070 

0071 

0072 

0073 

0074 


IV 0 Lfc/EL 21 GRNORN GATE * 7819* 14/47 

IF t I .EQ. J I X J1 ( 1 » J ) • ONE 

__IF J_i .IE. J) XJiCl.J) * ZERO _ . 

100 CONTINUE 

X J1C 2. 1 ) - -2 * RAW(l) 

XJ1C3.1I * -34R AM C 2 I ♦ 6*RAMI1)**2 
XJ1 C 3. 2 1 « -3 * R AM C 1 1 

X J 1 ( 4# 1) ■ -4*R AM ( 3 ) * 124RAWC2) *RAMC 1» - 12*R AM C 1 1 ** 3 
XJ1U.2) b*RjAM(ll*42 _ 

X Ji t 4. 31 - -4 * RAmCII 
C 

C INITI AL 1ZE 02 

C 

1)0 1C1 I * 1.4 
DO 1C1 J » 1.4 
IF (I .60. J ) GO TO 101 
XJ2CI.J) = ZERO _ 

101 CONTINUE 

XJ2C1.1I * ONE _ _ 

IF C XN . fcU . l.T GO TO 8 

7 XJ2C2.2) = 1. / CCENRLC2) » XN / CxN - 1.0) I 

XJ2C3.3) - ONE 

XJ2C4.4) * XJ2(2.2J**2 _/. 3.0 __ 

GO TO 9 

C _ _ 

C CALCULATE CHI -SQ JAR E TEST AND PARAMETERS 

C 

8 XN * XN ♦ 1. 

GO TO 7 

9 " CALL TRI PL E I XJ1. G.DJM) 

CALL TRI PLE( XJ2.DUM.SIGI ) 

CALL DMI NV (SIGI.4.DET.L.M) 

CALL RHAT2CM.SIGI.R) 

CALL AH4TC SlGl.R. A) 

CALL QUA T(XN.H.A.Q) 

CALL DGMPKDC R.H.THETA.4. 4. 1) 

TVAR « JEXP(THETaC2>) __ ___ 

WRITECb. 123) THETAC1 I.TVAR.O 

MR I TE ( 9. 124 ) THETA(l). TVAR. Q 

123 FORMAT ( /"/» T2 5. ' PARAMETERS : MJ * ' . E 1 5. 5. 10X. ' S I GMA *• . E 1 5 . 5 

£ / /» T39. * ***( CHI-SQUARE VALUE ) *** *»E15.5) 

124 FORMAT C / / . T2 5. ' NORMAL PARAMETERS: MU= '.Fb^^X,* SIGMA* *,Fb.2 

£ //,T39»* *♦*( CHI -SQUARE_ VALUE )*** '.F10.3) 

‘ RETURN 

END 
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: URTRAM IV 

G LEVEL 

21 GREXPQ 

DATE - 

78192 

0038 


CALL TRlPLEtXJl*G#SIGI I 




0039 


CALL DMlNVtSIGl»**D£T#L#N> 




00*0 


CALL RHAT1 ( W#SIGI» R) 




00*1 


CALL AHA T t S I Gl » R» A ) 




00*2 


CALL UH4TCXN#H, A, 0 1 




00*3 


CALL DGMPRDl R# H#THET A# 1 1 




00** 


XLAMJA « 1. / THETAI1I 




00*S 


WRIT E( 6» 123) XL AM CA» 0 




00*6 


WRITE! 9* 12* 1 XL AMDA# 0 




00*7 

i 23 

FORM AT( / /# T2 5» * PARAMETERS -* 

LAMOA « 

* » fci.5 . 

3 #/ 



6 lit T39» ' ***l CHI “S OUARE 

VALJc 

) *** • 

# £ 1 5 . 5 > 

00*8 

12* 

FORMAT ( V /# T25# * EXPONENTIAL PARAMETERS 

: LAMOA 

- *»Fo.2 »/ 



t / / » T39# ' *♦*( CHI-SUJARE 

VALJE 

) *** • 

*FL0.3) 

00*9 


RE TORN 




0050 


END 






Effect of Correlated Observations on Confidence 
Sets Based Upon Chi-Square Statistics 

Summary 

This paper investigates how the presence of correlation in a multivar- 
iate sample effects the confidence coefficients of confidence sets based 
upon chi-square statistics. 


I. Introduction 

Basu et. al. (1976) investigated the effect that simple equi correlation 
within a multivariate normal sample has upon confidence sets based upon chi- 
square statistics. They suggested that their results could provide a useful 
application in the area of pattern recognition using remotely sensed LANDSAT 
data. However, several recent investigations have demonstrated that the equi- 
correlated correlation structure is not an appropriate model in the Landsat 
application. In fact, Tubbs and Coberly (1978) demonstrated that the correl- 
ation struction in the LANDSAT data is similiar to observations obtained 
from a stationary autoregressive process. In this paper, I have investigated 
the effect that autocorrelated data have on confidence sets based upon chi- 
square statistics. 


II. Basic Concepts 

Let X^,...,X n denote a sample of n p-dimensional normal observations 
with mean u and common positive definite covariance matrix E. Suppose 


123 
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that X ■ [X, ,X_...,X ] T and that 
12 n 

E[ (X - E(x) ) (X - E(X) T ] = r S l (1) 

n 

where is a positive definite nxn matrix, A H B denotes the Kronecker 

product of matrices A and B, and E(«) denotes the expectation operator 

Note, if the sample X-...X is random then T =1, where I is an identity 

in n 

matrix. 

Now suppose that the sample X^**.X^ * s a realization from a discrete 

stationary time series {X } with continuous density function f„(*). If T 

t X n 

denotes the autocorrelation matrix for n lags. 

That is , 


~ (p^j) ^ >J ” 1*2. . . ,n 


"u ' corr( W- 


( 2 ) 


It is well known [Fuller (1972) ] that there exists an orthogonal matrix 
U such that 


where 


u*r U ~ 2n D v ( 3 ) 

n “ x 

= diag (d 1 ,d 2 ,.. . ,d n ) 

di . f x (0) 

d n ' V" 

PTIk 

^2k = d 2k+l "* f X ^ “rT” '* k = 1 » 2 » • • • >(n-l) /2. 
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and 


nW- 


2 ->s s -h 

1 cos(2II/n) 

0 sin(2H/n) 


cos ( 211 — ) 
n 

sin(2n — ) 

n 


1 

cos(^— ^ 2lT/n) . . . 

,n-l 
cos - 

2n — ) 


n 

n 

n 

0 

sin(— 2n/n) . . . 
n 


2 I,s r> 


By letting 

Z = U*X (5) 

it follows that 


E[ (Z - E(Z) ) (Z- E(Z) ) T ] = D x B E. 

Furthermore, it follows that 

lc - _ i n 

Z, ■ n 5 X; X = - E X, 
i n i 


( 6 ) 


(7) 


T 

vhsrs Z *s [Z-...Z ] . Bie distribvttion for Z. is 
in J 

Z 1 ^ NCn^y^E) 

Zj -V N(^,djE); J = 2 , 3 , . . . ,n 


( 8 ) 
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where the symbol 'v* means "is distributed as", 
since 


The expectation of Z 


J 


* zero 


Now let 


i 

E.(Z.) - E( E (cos(J-i 2Ilk/n) X. ) ) 
J k=0 n K 


(9) 


or 


E(V (sin(^- 2IIk/n) X. ) ) 

k=0 n k 


1—1 n-1 , . 

p( E cos(*“- 2llk/n) ) or p( E sin("— — 2Hk/n) ) 


k=0 


n 


k=0 


= 0 . 


Q-^p) = n(X - p) T E" 1 (X - p) 

n — T -1 — 

Q = E (X - X) 1 E ^(X, - X). 

j=l J J 


( 10 ) 


If T = I , it is well known that 


Q x (p) * x 2 (p) 

Q 2 * X 2 (n-l)p 


( 11 ) 


where 


X 2 (v) 


denotes a chi-square distribution with v degrees of freedom. 
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However, if T is given by (l) we have 
n 

Q^y) - n(X-y) T E -1 (X-y) 

- r 1 (nYA) 

- (Z^E^) ) T Z" 1 (Z 1 -E(2 1 ) ) 

- d 1 (Z 1 -E(Z 1 ) ) T (d^T 1 (Z 1 -E(Z 1 ) ). 

Hence 


Ql(u)/'di ^ x (P)' 


Now consider « E (X -X)^ I (X.-X) 

J-l J J 


tr E _1 [ E (X.-X) (X,~X) T ] 
J-l J J 


= tr E“ 1 [ E X.X T - riCX*] 
J-l J J 


However, since U is orthogonal (lU) becomes 


n 


Qo * tr E_1 t Z Z Z T - nXX A ] 
J-l J J 


n 


-i “ *p 

tr E 1 E Z Z 1 ] 
J-2 J J 


n 


E Z T E _1 Z. 
J-2 J J 


" Yj 


J*2 


(12) 


(13) 


(1U) 


(15) 




1?8 


for W « Z, T (d r) -1 Z . We knew that 
J J J J 

p degrees of freedom and that W^, 
i ^ J * 2,3,. . . ,n . 


W has a chi-square distribution with 

J 

are independent for each 


III. Confidence Set for Mean 

Let H denote the null hypothesis that X, ...X is a random sample from 
o in 

a p-dimensiona3 normal population with E(X) * y, cov(X) •= E, The statistic Q^, 
as given in equation (10) is used to define a confidence set for the unknown 
population mean u. That is, let 

I e = (y: Q x (y) < x e 2 (p)> ( 16 ) 

2 2 
where x (p) I s the 100 e percentage point of x (p). Thus since 
e 

2 

Q 1 'v x (p) whenever K q is true, we know that 

P[y e I £ | H q true] * e. (17) 

Let denote the alternative hypothesis that the sample satisfies «auation 

(l). If is true, then find the value a such that 


P[y c i c | h x 


true] 


a. 


(18) 
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From equation (13), ve know that a must satisfy the following relationship 


x/fp) ■ 


(19) 


IV. Confidence Interval for the Dispersion Scalar 

Let X^...X n denote a sample from a normal distribution with mean u 

2 

and covariance matrix o E, where E is a known positive definite matrix. Let 

H denote the hypothesis that the sample is random and H. denote the 
o l 

hypothesis that the sample satisfies equation (l). If H q is true, then 


Qg/a 2 * X 2 


p(n-l) 


where is given by equation (10). Hence the interval 


(20) 


°£ ° 2 £ VX E ,p(n-l) 


( 21 ) 


2 

is a 100 e confidence interval for o . However, to find the confidence 
2 

interval for a when H^ is true, it is necessary to determine the distri- 
bution of Qg. From equation (l?) we obtain 

• 2 n 

■ E d W (22) 

J*2 J J 

where W^, for J ■ 2,3,<««,n are distributed as independent chi-squares with 
P degrees of freedom. The distribution for (22) can be expressed in the 
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following series representation [c.f. Kotz, Johnson, end Boyd (1967)]. 

m 

P[Q,/o 2 < y] - E c. 6(v ♦ 2 k; y/e) (23) 

2 k-0 


where G(v+2k;y/0) denotes the cumulative probability density function for 

a central chi-square with degrees of freedom v+2k, and c^, 0 are known functions 

of the dj 's, for J a 2,3,...,n. Hence, whenever is true, the confidence 
2 

Interval for o in equation (21) is given by a where a is the value which 
satisfied the following relationship 


where 


a - Z c. G(n(n-1) + 2k; y e/ 6). 
k-0 K 


(24) 


y e- X c!p(n-1)* 


V. Examples 

Suppose that X^.-.X^ are a realisation from a stationary auto- 
regressive process of order one with parameter Then the spectral density 

function is 

f.(w) - = (25) 

A 2 n (1+4 -2^ cos w) 

Hence 

d 2k " < 14 + 2 - 2 + cos(2kIT/n) )~ l k-1,] n-1/2 

d x - ( l -^)" 2 


(26) 


13 ) 


The (lvalues which satisfy equation (19) are given in Table 1 for 
c - .99, 95. 


TABLE 1 

a- Values for AR(l) Process 


p\4> 

i 


2 


5 


10 


.0 

.1 

.2 

9900 

.9795 

.9606 

9500 

.9222 

.8830 

9900 

.9760 

.91*75 

9500 

.9116 

.8529 

9900 

.9661 

.911*5 

9500 

.8096 

.2856 

,9900 

.9570 

.8623 

,9500 

.861U 

.6952 


.3 

.1* 

.5 

.9285 

.8776 

.8021 

.8298 

.7603 

.6728 

.895 3 

.0O9U 

.6838 

.7695 

.6598 

.5270 

.8071 

.631*6 

.1*171* 

.6337 

. UL85 

.261*2 

.670U 

.1*055 

.1682 

.U6U8 

.2363 

.0 823 


. 4 * 
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Prom Table 1, we observe that a 95 % confidence elipse is a 65.98)8 
confidence elipse if the sample X^...X r is a bivariate sample from an auto- 
regressive process of order 1 with parameter <|> = . 1* 


TABLE 2 

a- Values for AR(l) Process 


N 

p\t> 

.0 

.1 

.2 

.3 

.1* 

.5 

.8 

13 

1 

.9500 

.9326 

• 8759 

.7901 

.6913 

.5938 

.3896 


2 

.91^3* 

.8817 

.7768 

.6317 

.1*822 

.3539 

.1518 


5 

.9lUU* 

.8211 

.5822 

.3365 

.1666 

.0751* 

.0089 

25 

1 

.91U3* 

.871*2 

.7577 

.5996 

.1*386 

.3020 

. .0902 


2 

1.0000* 

.8935 

.61*52 

.3869 

.1998 

.0927 

.0091 


5 

1.0000* 

.75U7 

.331*1* 

.0931* 

.0178 

.0026 

.0000 

51 

1 

1.0000* 

.8859 

.6223 

.3550 

.1702 

.0712 

.0036 


2 

1.0000* 

.7850 

.3872 

.1260 

.0286 

.0050 

.0000 


5 

1.0000* 

.51*60 

.0933 

.0056 

.0001 

.0000 

.0000 

101 

1 

1.0000* 

.7822 

. 3811 

.1209 

.0266 

,001*3 

.0000 


2 

1.0000* 

.6123 

.11*53 

.011*6 

.0007 

.0000 

.0000 


5 

1.0000* 

.2932 

.0080 

.0000 

.0000 

.0000 

.0000 


* the specified level e = .9500 

p 

From Table 2, a 99 % confidence interval for o is a 19.98)8 confidence 
based upon a bivariate sample of 25 observations from an AR(l) process with $ = 
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It is well known in applications using atmospheric observations that the 
data are non-random and in fact are highly correlated. Very little research 
has been done in the area of determining the effect that correlated samples 
have upon statistical inference. In this paper, I have investigated the effect 
thSt samples taken from a stationary autoregressive process have upon the 
confidence regions for the parameters of a normal distribution. Tables are 
included for the effect that sampling from an AR(l) process have upon these 
confidence regions. 
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GENERATION OF RANDOM VARIATES FROM SPECIFIED 


DISTRIBUTIONS 


Summary 

Due to the complexity of many of the existing statistical problems 
associated with atmospheric variables, computer simulations have proved to 
be a very informative technique. However, due to the various types of 
atmospheric data, thus the different type of statistical distributions one 
can no longer perform simulations based solely upon normal data. So in 
anticipating this problem, this paper presents the computer software for 
generating both random and correlated data for several specified distributions. 
A brief explantion of the procedure is given along with the program documen- 
tation. 


I. INTRODUCTION 

In order to obtain insight into some of the statistical problems 
with atmospheric data, it is necessary to be able to simulate some of the en- 
vironmental situations. However, since most of the data are non-normal it 
is necessary to generate data from various specified distributions (e.g. Gamma, 
Beta, Negative Binomial, etc.). The purpose of this paper is to document the 
procedures used in generating both correlated and uncorrelated observations. 

The uncorrelated procedures have been documented in Newmann and Odell (1971). 
The correlated procedures have been compiled from numerous sources, however, 
Johnson and Kotz (1972) provide the primary reference. In this paper, I have 
included only a brief description of the statistical distributions. For a 
more detailed discussion see Falls (1971). 
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II. UNCORRELATED VARIATES 

All of the procedures listed here are transformations 
of independent random variates from a uniform U(0,1) distribu- 
tion. The pseudo-random number generator used is a congruent ial 
generator (IBM SSP RANDU) whose choice was based solely upon 
convenience. However, some additional testing will be necessary 
to determine if the pseudo-random variates procedures are satis- 
factory for our purposes. 

Continuous Distributions 

o 

2.1 Univariate Normal Distribution N(u ,g ) 

The Box-Muller transformation [ 1] has been used. It 
can be summarized in the following result. 

Result: 2.1 If u and v are independently distributed 

U(0,1) then, 

x = (-2 In u)# cos 2 nv 
y = (-2 In v)^ sin 2 irv 

are independent random variates with the standardized normal 
distribution N(0,l). 

Thus if u^....Ujj is a sequence of independent U(0,1) 
one can generate a sequence x^....Xjj of independent N(0,1) 
using the above procedure. Also if a , o is a fixed known constant* 
then y^ = ox^ + p, i=l,2,...,n is a sequence of independent 

p 

normal with mean = p, variance = a • 
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2.2 Multivariate Normal Np(y,E) 

Let x^....Xp be a sequence of p independent normals 

m 

with mean 0 and variance 1, then x = (x,,....,x ) is said 

• r 

to be multivariate normal with mean 0 and covariance matrix 
I (pxp identity matrix). However, if x ^ N (0,1 ) then 

r r * 

y = Bx + p has a multivariate normal distribution with mean= u 

T 

and covariance matrix z , where I = BB . Prom x we can find y 
for any specified real positive definite symmetric matrix £ . 
This follows from the following result. 

Result: 2.2 Let z be a real p.d. symmetric matrix. Then 

there exists a lower triangular matrix B with positive eie- 

T 

ments on the main diagonal such that E = BB . This is often 
referred to as the Crout factorization of £ . 

2.3 Gamma Distribution r(x,k) 

Let u^....u^ be a sequence of k independent random 
variables each having a U(0,l) distribution. Then 

k 

x = - 1/ x In n u. (2) 

i=l 1 

is a gamma with parameters \ and k. Note the chi-square 
distribution with n degrees of freedom can be obtained by 
letting k=n/2 and x=#. Also, if n is odd then y = x+w^ is 
chi-square with d.f.=n if x -v r(k = n-)£, X =#) with w^N(0,l). 
The exponential distribution with parameter X can also be 
obtained by letting k=l in (2). 
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2 A Beta Distribution B(p.q) 

If x 2 ^ r(l,p) and X 2 'v r(l,q) are independent then 
y = x^ / (x^+X 2 ) has a Beta distribution with parameters 
p and q. 

Discrete Distributions 

If the distribution function F x is known then we can 
generate pseudo-random numbers by using the inverse function 
. However, this procedure can be simplified by letting 
x be the random variate from F x which satisfied the relation 
F (x-1) < u < F (x) where u is a random variate having a 

A " A 

U(0,1) distribution. This procedure could be used to generate 
Binomials, since the distribution function for the Binomial 
is easily obtained. Included is a discussion of some other 
discrete distributions which can be generated without knowledge 
of P x . 

2.5 Poisson Distribution P( a) 

If x^....x^ is a sequence of N independent exponentials 
with parameter X, then a non-negative integer k such that 
<_ 1 and 1 is distributed Poisson with parameter X , 

where 

k 

2.6 Negative Binomial Distribution NB(p,N) 



The negative binomial distribution can be generated 
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from a mixture of a Poisson and a Gamma distribution. That 
is, let Y be distributed as a Poisson with parameter e, where 
0 is a random variable from a Gamma distribution with parameters 
A,R. Then Y is distributed as a negative binomial with para- 
meters p = A/(l+ a) and N=R. 

Ill . CORRELATED VARIATES 

Continuous Distributions 


5.1 Correlated Multivariate Normal Distribution CNORM (u,E,A) 


Let Z Q ,Z^ r ..,Z^ be a sequence of N+l p-dimensd onal 
independent multivariate normals with common null mean vector 
0 and pxp covariance matrix z. Then 

X i = a i Z o + ^ 1 - a i^ Z i + p for i=l ,2 , . . . ,N 


are correlated multivariate normals with mean vector p and 
dispersion matrix A fii where 0 denotes the Kronecker product 
of A and z , that is 


A » Z = 
nxn pxp 


a-^E a^^ £ .... 

a2“^l' 1* . * . . 



(np x np) 


a ln Z 


a2nE 


a E 
nn 


and A is an N x 


N matrix where the i 



element of A is 


«i«d i ...,n 

1 i=3 


From the dispersion matrix A ® Z we have that 


i/d 

i-d 
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COV (x^Xj) = « i «. I 


d' " i d 

* E 


Hence the correlation matrix between vector X.pX^ is 


CORR(X.,X d ) = 


where 1^ is a pxp identity matrix, 
the univariate case. 


i/,1 

i-d 

When p is 1 we have 


5.2 Correlated Univariate Gamma Distribution r( X,R, A) 


Let Z Q , Z^, . . . .Z n denote a sequence of independent 
variables having the following Gamma distributions 

z 0 - r(»,8 0 ) 
z i »r(x,R i -R 0 ) 

Let X i =Z Q +Z^,i=l,2, . . . ,n, then X^....X n is a sequence of 
correlated Gamma variables where X^ r(A,R^) and the correlation 
between X^ and X^ is 


where a 


id 


CORR (X i ,Xj) = a id 

i.1. 

is the id element of the nxn matrix A and 


a 



1 if i=d 

R o # 

iqiq > if W 


( 
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3.3 Correlated Beta Distribution 8 (p.o t A ) 


Let Z Q Z^.... Z R be a sequence of independent chi- 
squares with degrees of freedom df=v^ (Gamma with A=l, 
^i =v i // 2^ for . . ,n. Let 

i=l ,2 , . . . ,n 


n 


X. = Z. / ( r Z.) 

1 1 d=0 3 

then the X^'s are correlated Beta with parameter (p^q^) 

n 

where p. = v . /?. and q. = p - p. where p * l p H . 
11 1 1 j>=o ' 3 

Then the correlation between X. and X. is given by 

J 

CORR (X.,X d ) = a id 

and 

a. 


id 



i"0 


7 p i p j ^ 

V (p-p^Cp-P-) 


Discrete Distributions 

3.4- Correlated Poisson P(X ,A) 

Let Z Q ,Z^, . . . . ,Z n be a sequence of independent Poisson 
with parameters ,i=0,l ,2 , . . ,n, then 

X i * Z o + z i 

is a sequence of correlated Poissons with X^ v P(a^) 

Ai = Ci + C Q , i=l ,2 , . . . ,n and the correlation between X^ 

Xj is given by 

Corr (XjXj) - a.j 
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and 



i/,1 


IV. CONCLUSIONS 

The purpose of this paper is to document the 
procedure used in programming uncorrelated or correlated 
number generators for various specified distributions. 

The results are fairly well known and should prove to be 
satisfactory for most simulation needs. As mentioned in 
the introduction, the procedures are dependent upon the 
choice of the pseudo-random number generator selected, and 
hence the objective of the situation to be simulated may 
dictate changes in the random number generator. A simple 
package is presented which would hopefully satisfy the 
needs of those researchers interested in generating numbers 
from the statistical distributions given. 
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APPENDIX A 

JOB CONTROL PARAMETERS 


CARD 

COL 

DESCRIPTION 

1 

1-5 

NREPS - Number of sets of numbers to be 
generated (15) 

(215) 

6-10 

IX - Seed for random number generator. 

(15) IX»0, then program will initiate 
using CPU clock 

** Note 

the following set of cards are repreated NREPS times 

2 

1-5 

NOB - Number of observations to be 
generated ( 15 ) 


6-10 

ITYPE - Type distribution to be generated (15) 

1 - Normal 4 - Poisson 

2 - Gamma 5 - Negative Binomial 

3 - Beta 6 - Binomial 


11 

ICOR = 1 correlated data (11) 

* 0 uncorrelated data ( 11 ) 


12 

ISTAT= 1 RJut Statistics (11) 

= 0 N • r it 


12-13 

IUNIT= 0 D ot output generated data ( 12 ) 

/ 0 Generated data output on 

external device # IUNIT 

*** Note the followini 
selected on Card j 

r cards depend upon the distribution 

f 2 . 

- NORMAL - 


3 

1-5 

NV * Number of variates (NV=2=bivariate 
normal ( 15 ) 


6-10 

KEY= 0 Standardised normal mean = 0 

variance = 1 



KEY =1 Read Mean, Variance ( 15 ) 


i 
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IF KEY - 1 Read following cards 
CARD COL DESCRIPTION 

4- (16F5.0) Y(I),I=1,.NV Mean vector 

5 (16F5.0) S(I) t I«l*NV**2 Covariance matrix 


** OF ICOR » 1 on card 2 read following for correlated case 


6 


Correlation factor (see page 4 ii) 


7 


Means (same grouping as correlation 
factors) only need when NV=1 


- GAMMA - 
3 1-5 

6-10 


R1 Shape parameter (F5.0) 

XLAMDA Scale parameter (F5.0) 


** IF ICOR - 1 Read following 


4 + 


Correlation factor (page iii) 


- BETA - 

3 1-5 R1 Beta parameter (F5.0) 

6-10 R2 Beta parameter (F5.0) 

** IF ICOR = 1 Read following 

4 1-5 VND Parameter for Z Q (see page 9 )(F5.0) 

5 + V(I), same format as correlation 

factors (page aii) 

- POISSON - 

3 1-5 XLAMDA Poisson Parameter (F5.0) 

** IF ICOR *= 1 Read following 

4 + Correlation factors (page iii) 


ii 
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- NEGATIVE BINOMIAL - 


CARD 

COL 


DESCRIPTION 

3 

1-5 

P 

parameter 

(*5*0) 


6-10 

N 

parameter 

(15) 


**IF ICOR * 1 Read following 
4 + Correlation factors (page iii) 

- BINOMIAL - 


No additional inputs needed. 


*** The following cards are used to define the A-matrix 
used in defining correlated observation 


- CORRELATION FACTORS - 


1 (13) 

2 (1615) 

3 (16*5.0) 
example 1 


example 2 


NG= Number of groups 1 1 NG 1 NOB 

NOL(I) .I=1,NG Length of each group 
N0L(1)+ N0L(2 )+...+ NOL(NG) « NOB 


VALUE (I),I=1,NG, A value for each 
group 


NOB-25 NG«1 N0L(l)-25 
VALUE(l)=.8 

the CORR(X. ,X,)«(.8)x(.8) - .64 

1 U 

CARD 

1 m 

2 m25 

, m. 8 


NOB =25 NG=2 

VALUE (1)=.5 
then CORR(X i ,X^ 


CARD 

i m 


> 


N0L(1)«10 N0L(2) 

VALUE (2)-. 8 

.25 i,d ilO 

< .40 i < 10, o' > 

.40 J < 10, i > 

^.64 i, j >10 


=15 

10 

10 


2 WV10WW15 


3 m. 5WV.8 


(Note : Y denote blank column) 
iii 
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MAIN 

SUPER 


BETA 

GAMMA 

BINOM 

NORMAL 

POISSN 

NEGBIN 

CBETA 

CGAMMA 

CNORML 

CPOISN 

CNEGBN 

PRINT 

STATS 

RANDU 


PROGRAM DESCRIPTION 

main program to read in job parameters 
supervisor routine to direct the generation 
of data, computation of statistics and printed 
output . 

generates independent Beta variates. 

" " Gamma variates. 

” ” Binomial variates. 

w M Normal variates. 

" " Poisson variates. 

" " Negative Binomial variates. 

” correlated Beta variates. 

,s Gamma variates. 

" " Normal variates. 

" " Poisson variates. 

" " Negative Binomial variates, 

prints generated values and output on specified 
unit. 

calculates statistic for generated values, 
generates random uniform variates. 


iv 
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SUBROUTINES NEEDED BY A GIVEN ROUTINE 


CNORML 


RANDA 

NORMAL 

CNEGBN 

- 

GAMMA 

POISSN 

CBETA 

- 

GAMMA 

C GAMMA 

- 

GAMMA 

CPOISN 

- 

POISSN 

NORMAL 

- 

RANDU 

BETA 

- 

GAMMA 

GAMMA 

- 

RANDU 

NORMAL 

BINOM 

- 

RANDU 

GMETRC 

- 

RANDU 

POISSN 

- 

GAMMA 

NEGBIN 

- 

GMETRC 

RANDA 

- 

RANDU 

RANDU 

- 


STATS 

- 


PRINT 

- 


MAIN 

- 

SUPER 

SUPER 


CNORML 

CBETA 

C GAMMA 

CNEGBN 

CPOISN 

NORMAL 

GAMMA 

BETA 

NEGBIN 

BINOM 

POISSN 

STATS 

PRINT 


v 



PRINT 


END OF 
STEPS 


GO BACK TO 
SUPERVISOR 


END 




















Arr 
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l)! M PK!OM X , Y( 100) ,S ( 1 00) ' /( 100) 

COMM0M/O/ Y.S,/ 

('OMMOM / A / IX .‘IV. 91 ,XLAM0A,92«P «M 

I X='bl 

RFAnls.inn) m*>fps 

KO 99 11=1 .MPFPS 

PF'AD<*>. 100) MO. ITYPF. IUmIT 

CALL TYPF(X.MO. ITYPF ) 

CAL L OP I M T ( X . MO . I DM I T * I TY PF • 1 1 ) 

QQ COMTIMHF 
inn F OP MAT (IIS) 

STOP 
f NO 


SiiPPOllT T KI F »FTA(X.MO) 
n iMFN<;inn x ( ? on ) • y ( i n o ) .s u on ) . ? < ioo) 

COMMOM/A/ IX, MV. PI , X L A M O A , P 2 » P » M 
covmom/p/ Y.S, 2 
XL A Ml) as 1 . 

CALL OAMMA(X.MO) 

P=»l 
Pl = -V 

CALL GAMMA <Y. MO) 

P l = p 

oo 1 T = 1 , MO 
XX=X(I)/(X(I)+Y(I>) 

1 X C T ) = X X 
RtTUPM 
EMO 


SOPWOmiMF GAMMA ( X, MO) 

dimension x(Mn).Y()no) 

COMMOM/P/ Y.S./ 

COMMOM/A/TX »m V» 01 ,XLAMOA,P?,P.M 
XL=-1 .*< 1 ,/XLAMDA) 

K=P1+.S 

u — *) 1 _l f 

DO 99 T 1=1 , MO 
X X=1 , 

0 (( 1 1 = 1 . K 
I M= I X— ( FX/10)*10+1 
DO 2 J=).IM 

■? CALL PAMOU ( IX .IX.YFL) 

1 XX=XX*YFL 
99 X (I I) =XL*ALOG (XX) 

IF(W.LF.O) PPTLtPM 
NV = 1 

CALL mopmaL (Y ,MO. MO .0 ) 

DO 9M T=] ,Nn 
9B X (I )=X( I) +Y (T )*Y< l) 

R FT URN 
EMO 

ORIGINAL PAGE Id 
OF POOR QUALITY 

SIRWOmtimE pIMOM(x.mO) 

DIMENSION X(M0) 

COMMOM/A/ I X»MV .PI ♦ XLAMOA ,R2,P.N 
DO 1 1=1, MO 
lM=IX-( IX/] 0) *1 0+1 
DO 2 K=1 . IM 

? CALL PAmO(J( IX . I X. YFl) 

1 X(I)=YF1 

pf.tupm 

E Ml) 
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AT'T’l' !)!.'•: 


S • m po i it i * •(■* mo>vmai ( x ,wn. mvukf V) 

n i mfms t * 1*' * ( ^on ) . y ( i oo ) . s ( i oo ) « / ( i no > 

Common/ f’ / y , s . / 

COMMON/ A/ IX .MV. ‘M .XL AMD A.*? 

* . P.N 

NVV=NV < »MV 
P ?—h . POO 1 MS 

C OF NFP AT F mv'O TMOFPFMOF'MT IIMIFOPM (0.1) 

DO P 1 = 1* mV n 
I m= lx - ( I < / 1 0 ) * l ('♦ 1 

DO ( k' = l,I..| 

0 CALL PAMI1!) ( I X , I X . YF|.) 

? X ( I )=Yf| 

C T P A M S F 0 P 1 TO MVO MOWMAL (0,1) 

J =M VO /P • 

DO '» T ” 1 . .J 

A=S'JP T ( -P .* A 1.00 (X (?* 1-1))) 

x ( ?* I -i ) = A*cno(P,>ox ( i»?) ) 

4 > (?*I ) = A*SI M(‘>P*X ( I »? ) ) 

IF(KFY.FQ.O) PFTllPM 
FPI TF (ft . 00 ) 

POO F 00 M A T ( • MF A* AMO COV AM I A*.|CF« » // ) 

PFAo<S. 100) ( Y( T) .1=1 .MV) 

WWITL(«..?01) (Y (T ) . 1 = 1 .MV) 

I ' F A 1 1 { S , 100) (S ( I ) . 1 = 1 . MV V) 

VDMTh (ft.POl ) ( S < I) • T = 1 . M VV ) 

?01 F OP MA r ( 1 0 X » 1 (' F 1 0 . 1 ) 

IF (NV .OT. I ) 00 TO ft 
DO S t=i*mo 

5 X ( I ) =S ( 1 ) »X ( I ) +Y ( 1 ) 

RFTUHm 

ft M=MV 

DO 11 J = 1 .MV 
f • o 1 1 1 = 1.0 
11 l ( ( j- 1 ) *M+I) = 0. 

/ ( 1 )=SOPT (S (1 ) ) 

001 3 1 = P * N 

n i ( i ) =s ( t ) /z ( i ) 
noio ,j=p.m 
DO l v t=i.m 
SI)M=0 . 

IF(I-.J) IP. IS, 17 
IS M=.)-i 

DO 1ft *=1 .M 

lft S(JM=5 i)M 4'/ ( ( K-l ) «M + T )*«? 

I ( ( J— 1 ) *m+i )=SOPT (S ( ( J-l ) *N + T )-SIIM) 

GO TO 19 
17 M=J— 1 

noift k = i . m 

iq SlJM = Si 1*1+ ( i ( (K-l )*N+T) *? < (K-l ) *M+J) ) 

/ ( ( J-l ) <>M + I )= (S ( ( J- 1) *M + I )-SUM) // ( ( J- 1 ) »M+J ) 

1 9 CONTI Ml (F 

DO 7 T = 1 * MO 
DO 7 K=1 ,NV 
SMI-1) «MV + K) =Y (* ) 

DO 7 ,J=1 .K 

7 S ( ( 1-1 ) *MV+K) =S ( ( T-l ) *NV+K) + 7 ( ( J-l ) *MV+K) 

* *X (MVM- (NV-J) ) 

DO H 1=1, MO 
do y ,)=1.MV 

Q X ( ( J-l ) *M0+ I ) =S ( ( 1-1 ) *MV+J) 

PFTUPM 

inn FOHMAT(l^Fh.O) 

E NO 
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Apr 1 !' in.'; 


SUuwnnTTMr 0*4' TP C(XX) 
f OMM()m/ A/ I <• MV *P1 « XLAM f )A .P£, P.M 
X X = 0. 

P = . b 
OMf-=l . 
o=nNfc*-p 
no i T-i.M 
CAI.L OAMOIU JX .1 X • t n 
shm=o. 

J=l) 

(JO=P 
T J=J + 1 
oo=oo*n 
SllM=S'|M+no 
IF (bUM-11) ?. 1 . 1 
7 IF ( J.I.T.IO) r,0 TO 3 
\ XX=XX+.J 

pftuwm 

F Mil 


Si 10 00 1 IT IMF PO IS^M (X ,f.JO) 
UI m FMSTOM X(M(>) 

COMM0M/ a / IX.MV.P] . XLAM[)A.P?«P«M 

msum=o 

kl=l. 

1 SllM = U . 

J = 0 

? CALL. OAMM A ( xx . 1 ) 

J=J + 1 

SlH = SllM + X X 

I F ( SIIM.I.F . 1 .) 00 TO ? 
x IMStlM) =.l-l 

NSUH = MSIIM+1 

IF (NSlIM.l F.MO) 00 TO ] 

kf TIIRm 

FMD 


FHpRO'ITTmF TYPF (X .mo, ITYPF) 
DIMtNSTOM X (MO) 

COMMOM/A/ I X . N|V . P 1 « XLAMOA .R2.P.M 
(■>() TO (l.?.T.4»b*0) » I TYPF 
1 PFAO (S • 1 OO ) MV.KFY 
M VV=NV*MV 
MV/0=NV*M0 

CALL MOPMAL (X .MO.MVO.KFY) 

PPTUPM 

? PFAr)(S. 1 1» 1 ) Pl.XLAMOA 
CAU. oamma(x.mo) 

WFTUHm 

1 PFAlHS.ini ) R I ,P? 

CALL qFTA(X.MO) 

PF TUPM 

A PFAtMS, 1 01 ) XLAMOA 
CALL ROTSSO(X.MO) 
k FT (HIM 

S CALL WTMnM(X.NO) 

Hf.TUPn 

X, RFAD(S,1PP) P.M 
CALL MFGPIN(X.MO) 

PFTIIPO 

1 Of) forma t ( ?\ S> 

1 n i k ORM A T ( ?F S. 0) 

5 OP F OPMA T ( F 1 ^ ,0 . 1 b> 

F rjf) 


OF p- ;n 


"'Atsfe 
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Am: ;> t 


4,lin KOI IT T ►' F Mf- ,,n l N ( X ,M()) 
f 1 1 4' (• M s I ON X (MO ) 

C OMMOM/ A/ I x . mV «R 1 «XL AMpA ,P ?» P, N 

pn 1 (= 1 . mo 

CAL L GMF T«C( XX) 

x ( I )=xx 

RETURN 

FNR 


SUBROUTINE WAVJOIM IX, I Y. YFL) 

I Y= 

If ( IY)S.*.A 
s IY-IY«.? 14 74 h;^4 7+1 
A Y K I . = I Y 

YFL=YFL*.4*FvSM 3E-9 
RETURN 
F MO 


SI'OROHT IMF PP I M T (X ,MO, JIJ, I T . I I ) 

DIMENSION X(NO) 

COMMON/ />/ I X .MV, PI .XLAMf'A.pp.P.M 
N VO=MV tf F'O 

00 TO ( 1 .P«^.4»o,f,),iT 
1 WRITK (0.1 DO) I I. MV 

WRITE (o.? »0) (X (I ) . 1=1 . NVO) 

IF( IU.fO.?) WRITE (9,300) (X( I) ,1=1 .NVO) 

RETURN 

P WPITE(O.lOl) II ,P 1 • XL AM OA 
WRITE (S.POO) X 

IF(HJ.F0.2) WRITE (9,300) X 
RFTURm 

3 WRITE (0,1 0?) II ,R1.R? 

WRITE (A.POO)X 

IFIIU.F0.2) WRITE (9,300) X 
RETURN 

4 WRI TE ( 0 ,1 cm II .xlamoa 

WRITE (0. POO) X 

lF(IU.FO.P) WRITE (9,300) X 
RFTURm 

5 WRITE ( A » 1 04) II 
WRITE (A.POO)X 

IEUU.FO.?) WR I TF (9,300 ) X 

return 

A WRI TF (A.l O'?) II , R , M 
W°I TF (A. POO) X 

inilJ.FO.(») WRI IF (9,300) X 
RETURN 

1 no FORMAT(//,i normal DATA FOR RF P .= • » Is ,//, • NO. OF V AR I A TE S= * . I A . 

* //) 

ini FORMAT(//.» GAMMA DATA FOR REP. = • , I 6, // . • N=»,F10.3» 

* • L AMO A = » .F) 0.3,//) 

10? FORMA T (//, • P FT A 04TA FOR RFP. =» , I A, // , i ALPUA=« ,F ] 0 .2 , * RET A= • , 

* P10.P.//) 

103 FORMAT!//.* POISSN DATA FOR REP , = • • I A • / /. * LAMD A= t . F 1 0. ?. // ) 

104 F ORMA T ( / / * * UNIFORM (0,1) RATA FOR INTEGRAL TRANSFORM PEP=*,IA,//) 
JOS FORMAT (//.* MEG. RlNOMIAL OAT A FOR RE p » =• » 16 ♦ // » 

* i P=i,F10.3,t M=»,I10,//) 

POO F 0«MA T ( 1 0 X* 1 OF 1 0. 3) 

300 FMRMAT(G0F10,3) 

F NO 


